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Abstract 

An  upward  Vertical  Flow  Treatment  Wetland  (uVFTW)  has  been  designed  to  use 
anaerobic  and  aerobic  microbial  processes  to  bioremediate  groundwater  contaminated 
with  chlorinated  hydrocarbons.  Hydraulic  short-circuiting  has  been  a  problem  with 
uVFTWs.  Corbin  (2008)  estimated  that  for  a  uVFTW  constructed  at  Wright  Patterson 
AFB  to  treat  contaminated  groundwater,  groundwater  flowed  through  less  than  50%  of 
the  wetland’s  volume,  and  that  the  actual  mean  residence  time  (1.38  days)  was 
significantly  less  than  the  8.75  days  that  would  be  achievable  with  uniform  flow  through 
the  wetland  cell. 

The  objective  of  this  research  is  to  investigate  how  the  hydraulics  of  uVFTWs 
affects  treatment  efficiency.  A  sub-objective  is  to  propose  uVFTW  design  strategies  that 
can  be  used  to  maximize  treatment  efficiency.  To  accomplish  this,  a  model  of  a  three- 
layer  uVFTW  that  couples  hydraulics  and  degradation  kinetics  was  built  using  the 
MODFFOW -based  Groundwater  Modeling  Software.  The  model  was  applied,  using  data 
collected  by  Corbin  (2008)  and  literature  values  for  degradation  rate  constants,  to 
estimate  the  effectiveness  of  various  engineering  solutions  aimed  at  improving  overall 
treatment  efficiency. 

The  results  indicate  that,  compared  to  a  baseline  model,  which  was  constructed  to 
approximate  the  existing  Wright  Patterson  AFB  uVFTW  and  had  a  mean  residence  time 
(tau)  of  6.19  days,  all  of  the  proposed  engineering  solutions  would  be  effective  in 
increasing  the  mean  residence  time.  Based  on  the  actual  flow  through  the  wetland,  and 
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assuming  the  entire  wetland  volume  was  used  (i.e.,  no  hydraulic  short-circuiting)  the 
nominal  hydraulic  residence  time  of  the  wetland  cell  would  be  9.50  days.  Using  model 
simulations  to  assess  the  impact  of  various  designs  aimed  at  improving  hydraulic 
performance,  the  mean  residence  time  increased  10%  (tau  =  6.89  days)  to  180%  (tau  = 
17.5  days)  from  the  base  line  tau  of  6.19  days. 

The  effect  on  treatment  efficiency,  however,  was  less  clear.  While  there  was  a 
significant  improvement  in  treatment  efficiency  when  comparing  the  proposed 
engineering  solutions  and  the  baseline  model,  almost  the  entire  improvement  occurred 
with  a  very  small  increase  in  mean  hydraulic  residence  time,  from  6.19  days  to  6.89  days. 
Treatment  efficiency  was  essentially  unchanged  as  mean  hydraulic  residence  time 
increased  from  6.89  days  through  17.45  days.  This  is  due  to  the  relation  between  the 
degradation  time  constant,  which  is  determined  by  the  degradation  rate,  and  the  hydraulic 
residence  time,  tau.  For  the  baseline  conditions,  the  first-order  degradation  rate  constant 
for  the  influent  contaminant,  tetrachloroethylene  (PCE),  was  0.4  d’1.  This  corresponds  to 
a  half-life  (i.e.,  the  degradation  time  constant)  of  PCE  of  1.7  days.  Since  the  degradation 
time  constant  is  comparable,  but  smaller  than,  the  baseline  tau  of  6.19  days,  an  increase 
in  tau  leads  to  improved  destruction  efficiency.  However,  when  tau  is  sufficiently  large 
compared  to  the  degradation  half-life,  increases  in  tau  have  negligible  impact  on 
efficiency. 

This  study  demonstrates  that  a  quantitative  understanding  of  contaminant 
degradation  kinetics  is  important  in  making  decisions  regarding  constructed  uVFTW 
hydraulics.  Hopefully,  this  study  will  serve  to  guide  designs  of  future  uVFTWs. 
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MODELING  VERTICAL  FLOW  TREATMENT  WETLAND  HYDRAULICS  TO 


OPTIMIZE  TREATMENT  EFFICIENCY 


I.  Introduction 


1.1  Background 

An  upward  Vertical  Flow  Treatment  Wetland  (uVFTW)  has  been  designed  to  use 
anaerobic  and  aerobic  microbial  processes  to  bioremediate  groundwater  contaminated 
with  chlorinated  hydrocarbons  (Amon  et  al.,  2007).  In  the  design  that  has  been  used  at 
Wright  Patterson  AFB  (Amon  et  al.,  2007)  the  tetrachloroethylene-  (PCE-)  contaminated 
groundwater  is  pumped  into  the  bottom  gravel  layer  of  the  wetland  cell.  The  water  then 
flows  upward  through  three  different  cell  layers,  where  bacteria  biodegrade  the 
contaminants.  Treated  water  then  exits  over  a  weir  at  one  end  of  the  rectangular  cell. 

The  design  objective  is  to  maximize  treatment  efficiency. 

Hydraulic  short-circuiting  has  been  a  problem  with  uVFTWs.  Corbin  (2008) 
estimated  that  for  the  uVFTW  at  Wright  Patterson  AFB,  groundwater  flowed  through  less 
than  50%  of  the  wetland’s  volume.  Corbin  (2008)  also  calculated  that  the  actual  mean 
residence  time  (1.38  days)  was  significantly  less  than  the  8.75  days  that  would  be 
achievable  with  uniform  flow  through  the  cell. 

The  understanding  of  the  hydraulics  of  constructed  wetlands  and  the  processes 
whereby  contaminants  are  degraded  has  improved  dramatically  in  the  last  decade.  Initial 
models  of  constructed  wetlands  were  little  more  than  “black  boxes,”  giving  little 
understanding  of  the  biological  and  chemical  processes  occurring  (Langergraber  et  al., 
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2009).  Recent  mechanistic  models,  although  quite  complex,  are  capable  of  simulating 


multiple  biological  and  chemical  processes  (Langergraber  et  al.,  2009).  Few  of  the 
numerical  models  currently  in  use  deal  with  subsurface  flow  in  constructed  wetlands,  and 
most  of  those  simulate  horizontal  flow,  rather  than  vertical  flow  (Langergraber,  2008). 
Upward  vertical  flow  treatment  wetlands  are  an  even  smaller  subset  of  vertical  flow 
wetlands,  and  there  have  been  very  few  studies  modeling  their  hydraulics. 

1.2  Research  Objective 

The  objective  of  this  research  is  to  investigate  how  the  hydraulics  of  uVFTWs 
affects  treatment  efficiency.  A  sub-objective  is  to  propose  uVFTW  design  strategies  that 
can  be  used  to  maximize  treatment  efficiency. 


1.3  Specific  Research  Question 

1)  What  design  parameters  have  the  greatest  impact  on  treatment  efficiency?  Would 
altering  the  ratio  of  cell  length  to  width  to  depth,  relative  layer  thickness,  or  influent 
piping  configuration  improve  treatment  efficiency? 

2)  Based  on  the  answer  to  question  1,  what  design  changes  can  be  implemented  to 
maximize  treatment  efficiency?  How  much  of  an  effect  on  treatment  efficiency  would 
these  design  changes  have? 
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1.4  Research  Approach 

1)  Conduct  a  literature  review.  In  addition  to  seeking  information  on  whether  this 
has  been  tried  or  modeled  before  (i.e.,  how  have  models  been  applied  to  predict  treatment 
efficiency  as  a  function  of  hydraulics  and  degradation  kinetics  in  constructed  wetlands?), 
the  literature  review  will  be  used  to  determine  typical  values  for  environmental  and 
design  parameters,  as  well  as  typical  treatment  efficiencies.  The  literature  review  will 
also  include  case  studies,  where  the  impact  of  hydraulics  and  degradation  kinetics  of 
constructed  wetlands  (and  particularly  uVFTWs)  on  overall  treatment  efficiency  has  been 
assessed. 

2)  Build  a  generic  model  of  a  three-layer  uVFTW  that  couples  hydraulics  and 
degradation  kinetics.  For  simplicity,  biodegradation  will  be  described  using  first-order 
kinetics,  although  the  first-order  rate  constant  can  be  different  in  the  different  layers.  The 
validity  of  the  first-order  model  will  be  assessed  in  the  literature  review. 

3)  Use  the  model  to  conduct  sensitivity  analyses  to  see  how  varying  design 
parameters  (location  of  influent  piping,  baffle  placement,  cell  dimensions,  etc.)  as  well  as 
environmental  parameters  (degradation  rates,  existence  of  high  conductivity  flow  paths 
that  may  cause  short  circuiting),  affect  treatment  efficiency. 

4)  Apply  the  model,  using  data  collected  by  Corbin  (2008)  and  literature  values  for 
degradation  rate  constants  for  the  uVFTW  located  at  Wright-Patterson  AFB,  to  estimate 
the  effectiveness  of  various  engineering  solutions  at  improving  overall  treatment 
efficiency. 
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1.5  Scope  and  Limitations  of  Research 

The  most  significant  limitation  of  this  research  is  that  it  is  a  modeling  study,  and 
therefore  depends  on  the  applicability  of  a  number  of  simplifying  assumptions,  and  that 
the  only  data  that  are  available  come  from  the  single  uVFTW  being  used  to  treat 
chlorinated  hydrocarbons  in  existence.  Thus,  model  validation  may  not  be  achievable  in 
this  study.  However,  it  is  hoped  that  development  and  application  of  the  model  will 
provide  insights  into  how  uVFTWs  can  be  designed  to  more  efficiently  treat  water 
contaminated  by  chlorinated  hydrocarbons. 


1.6  Definition  of  Terms 

Term  Definition 

FWSW  Free  Water  Surface  Wetland 

HSFW  Horizontal  Subsurface  Flow  Wetland 

VFTW  Vertical  Flow  Treatment  Wetland.  The  suffix  “u”  or  “d”  will  be  used 

to  distinguish  between  upward  flow  (uVFTW)  and  downward  flow 
(dVFTW)  treatment  wetlands. 

CVOC  Chlorinated  volatile  organic  compound.  In  general,  they  have  high 

vapor  pressures,  low-to-medium  water  solubility,  and  low  molecular 
weights.  CVOCs  are  ground-water  contaminants  of  concern  because 
of  very  large  environmental  releases,  human  toxicity,  and  a  tendency 
for  some  compounds  to  persist  in  and  migrate  with  groundwater  to 
drinking  water  supply  wells  (U.S.G.S.,  2010). 

Halorespiration  Process  of  anaerobic  microbial  respiration  where  microbes  gain 
energy  through  the  use  of  a  halogenated  compound  (compound 
containing  fluorine,  chlorine,  bromine,  or  iodine)  as  a  terminal 
electron  acceptor  (Lorah  et  al.,  2007). 
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Term 

Redox 


Methanogenesis 


Porosity,  n 


Retardation 

Factor 


Critical  Flow 


Definition 

Oxidation-Reduction.  A  general  term  describing  reactions  involving 
the  transfer  of  electrons  from  the  oxidized  compound 
(reductant/electron  donor)  to  the  reduced  compound  (oxidant/electron 
acceptor). 

An  oxidation-reduction  reaction  involving  the  reduction  of  carbon 
dioxide  to  methane.  Under  sufficiently  reducing  conditions, 
methanogenesis  is  the  dominant  oxidation-reduction  reaction. 

Total  porosity  is  defined  as  the  ratio  of  the  volume  of  void  spaces  (Vv) 
to  the  total  volume  ( Vr)  in  a  porous  medium: 


Effective  porosity  is  the  percentage  of  interconnected  pore  spaces,  and 
is  always  less  than  total  porosity  (Domenico  and  Schwartz,  1990). 

Ratio  of  the  average  flow  velocity  of  groundwater  due  to  advection  to 
the  flow  velocity  of  a  sorbing  contaminant  (Domenico  and  Schwartz, 
1990): 


R  =  l  +  p-±Kd  (2) 

where, 

Pb  Soil  bulk  density  [M  L3] 

n  Porosity 

3  1 

Kd  Sorption  distribution  coefficient  [L  M"  ] 

Use  of  a  retardation  factor  implicitly  assumes  linear,  equilibrium 
partitioning  of  the  contaminant  between  the  sorbed  and  dissolved 
phases,  characterized  by  the  distribution  coefficient,  which  is  the  ratio 
of  equilibrium  contaminant  concentrations  in  the  sorbed  and  dissolved 
phases. 

In  a  uVFTW,  critical  flow  is  the  maximum  flow  that  can  be  applied 
per  acre  prior  to  fluidizing  the  bed,  and  has  units  of  [L3  T  1  L  2] 
(Pardue,  2005). 
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Term 

Hydraulic 

Conductivity 


RTDF 


Mean  residence 
time,  x  or  tau 


Variance 

(second 

moment) 


Definition 

The  capacity  of  a  medium  to  transmit  water.  Typically  determined 
experimentally  using  a  number  of  techniques  and  applying  Darcy’s 
Law  (Domenico  and  Schwartz,  1990): 


Q  =  KAVh  (3) 

where, 

Q  Volumetric  flow  [L3  T 1  ] 

A  Cross  sectional  area  perpendicular  to  the  flow 

direction  [L2] 

V/i  Hydraulic  gradient  [-] 

K  Hydraulic  conductivity  [L  T1] 

The  Residence  Time  Distribution  Function,  f(t),  describes  the 
probability  that  water  or  a  conservative  tracer  will  spend  a  given 
amount  of  time  in  a  reactor  bed.  Typically  displayed  graphically  with 
RTDF  on  the  y  axis  (units  of  inverse  time),  and  time  on  the  x  axis 
(Clark,  2009). 


The  mean  residence  time,  or  x  (tau),  is  the  first  moment  of  the  RTDF, 
f(t),  and  is  given  by  (Clark,  2009): 

T  =  /0°°/(0  x  tdt  (4) 


For  discrete  data,  this  equation  becomes: 


\'lma : 
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(5) 


The  variance,  also  known  as  the  second  moment  about  the  mean  of  the 
RTDF ,f(t),  describes  the  spread  of  the  RTDF  (Clark,  2009): 

0-2  =  /0°°(  t-T)2f(t)dt  (6) 


For  discrete  data: 
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Term 

Dimensionless 

variance 


tbar 


STELLA 


Definition 

The  dimensionless  variance  allows  an  evaluation  of  how  close  a 
residence  time  distribution  is  to  a  perfect  mixer  or  plug-flow  reactor: 


Lor  the  perfect  mixer  —  1,  while  for  the  plug-flow  reactor  <J?  =  0. 
Lor  most  reactors,  o%  is  bounded  between  0  and  1  (Clark,  2009). 

Nominal,  or  theoretical,  mean  residence  time  [T] .  Calculated  by: 

Vn 

tbar  —  "jj"  (9) 

where, 

V  Volume  of  treatment  wetland  [L  ] 

n  Porosity  [-] 

Q  Volumetric  flow  rate  [L3  T"1] 

Structural  Thinking  Experiential  Learning  Laboratory  with  Animation. 
Software  tool  designed  to  model  complex  dynamic  systems  (isee 
systems,  2010). 
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II.  Literature  Review 


2.1  Review  of  Treatment  Wetland  Designs  and  Purposes 

Treatment  wetlands  have  been  used  for  decades  to  treat  a  variety  of  wastewaters, 
including  municipal  wastewater,  wastewater  from  mining  operations,  industrial 
wastewater,  and  water  contaminated  by  livestock  operations  (Cole,  1998).  There  are 
several  different  types  of  treatment  wetlands,  each  having  advantages  and  disadvantages 
depending  on  the  type  of  wastewater  to  be  treated.  Treatment  wetlands  may  be  divided 
into  two  broad  categories:  surface  flow  and  subsurface  flow. 

In  a  Free  Water  Surface  (FWS)  treatment  wetland,  the  water  is  exposed  to  the 
atmosphere  as  it  flows  horizontally  through  the  cell  (Figure  1).  Contaminants  are  treated 
or  removed  through  sedimentation,  filtration,  interaction  with  the  aquatic  plant  root  zone, 
or  other  processes.  FWS  treatment  wetlands  are  typically  used  for  storm  water  runoff  or 
for  tertiary  treatment  of  wastewater.  Their  main  advantage  is  their  ability  to  handle  pulse 
flows  and  changing  water  levels.  Their  main  disadvantage  is  reduced  efficiency  during 
sub-freezing  weather.  An  unintended  benefit  is  the  creation  of  wildlife  habitat  (Kadlec 
and  Wallace,  2009). 

FWS  treatment  wetlands  have  been  used  to  treat  chlorinated  hydrocarbons 
(Kadlec  and  Wallace,  2009).  At  a  site  in  Minnesota,  contaminated  groundwater  flowed 
into  a  dredged  channel  which  conveyed  the  water  to  a  lake.  The  channel  was  filled  and 
planted,  converting  it  to  a  treatment  wetland.  At  the  Schilling  Farm  Project  in  Michigan, 
a  FWS  treatment  wetland  was  constructed  in  the  path  of  a  TCE-contaminated 
groundwater  plume  (Kadlec  and  Wallace,  2009). 
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Keefe  et  al.  (2004)  evaluated  the  fate  of  volatile  organic  compounds  in  a  wetland 
constructed  for  the  treatment  of  wastewater  treatment  plant  effluent.  The  wetland  is  a 
variation  on  the  FWS  system,  with  flow  proceeding  alternately  through  5  deep  zones  and 
6  shallow  zones.  Removal  efficiencies  for  the  compounds  studied  ranged  from  63%  to 
87%.  The  model  used  was  not  a  mechanistic  model,  so  the  specific  pathway  by  which 
mass  was  removed  was  not  simulated  (though  the  investigators  assumed  volatilization 
was  the  dominant  removal  pathway)  (Keefe  et  al.,  2004). 


Figure  1:  Typical  Free  Water  Surface  Treatment  Wetland 


In  a  subsurface  flow  treatment  wetland,  the  water  flows  entirely  beneath  the 
surface  and  all  treatment  occurs  beneath  the  surface,  as  well.  Since  the  water  is  not 
exposed  to  the  surface,  there  is  no  risk  of  human  contact,  and  more  polluted  influent  may 
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therefore  be  treated.  These  wetlands  are  commonly  used  for  secondary  treatment  of 
residential  wastewater,  either  for  individual  homes  or  small  communities.  Their  main 
advantage  is  they  can  operate  in  freezing  conditions.  Their  main  disadvantages  are  a 
propensity  for  clogging  of  the  subsurface  media  and  smaller  inlet  flow  rate  (Kadlec  and 
Wallace,  2009). 

Subsurface  flow  treatment  wetlands  may  be  further  divided  into  horizontal  and 
vertical  flow.  A  horizontal  subsurface  flow  treatment  wetland  (HSFW)  is  built  with  a 
sloped  bottom  to  ensure  flow  in  the  desired  direction  (Figure  2).  Influent  is  piped  into  the 
cell  just  below  the  surface,  at  the  up-gradient  end.  The  contaminated  water  flows 
horizontally  through  a  gravel  bed,  where  it  is  treated  through  contact  with  plant  roots. 

The  effluent  is  collected  near  the  bottom  of  the  down-gradient  end  (Kadlec  and  Wallace, 
2009). 
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Figure  2:  Typical  Horizontal  Subsurface  Flow  Treatment  Wetland 


A  vertical  flow  treatment  wetland  is  typically  understood  to  mean  downward 
vertical  flow  (dVFTW).  In  this  type  of  system,  the  treatment  wetland  is  “pulse  loaded” 
with  a  large  influx  of  contaminated  water,  which  is  then  treated  as  it  moves  downward 
through  the  cell  (Figure  3).  This  type  of  system  provides  a  higher  level  of  oxygen  than  a 
horizontal  subsurface  flow  system,  and  is  therefore  capable  of  oxidizing  ammonia. 
Downward  VFTWs  are  capable  of  handling  wastewater  containing  high  levels  of 
ammonia,  such  as  landfill  leachate  and  food  processing  wastes  (Kadlec  and  Wallace, 
2009). 
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Figure  3:  Typical  Downward  Vertical  Flow  Treatment  Wetland 


Upward  VFTWs,  such  as  the  one  constructed  at  Wright  Patterson  AFB  in  2000, 
are  a  more  recent  design.  The  impetus  for  this  specific  VFTW  was  the  work  of  Lorah 
and  Olsen  (1999),  who  observed  that  trichloroethylene  (TCE)  and  1, 1,2,2- 
tetrachloroethane  (TeCA)  in  a  plume  of  contaminated  groundwater  in  a  natural  wetland 
underwent  anaerobic  reductive  dechlorination  to  the  daughter  products  1,2- 
dichloroethylene  (cDCE),  vinyl  chloride  (VC),  1,1,2-trichloroethane  (1,1,2-TCA),  and 
1,2  dichloroethane  (1,2-DCA)  as  the  plume  proceeded  along  an  upward  flow  path  from 
an  aerobic  sand  aquifer  through  an  anaerobic  region  of  a  wetland  (Lorah  and  Olsen, 
1999).  Conditions  in  the  anaerobic  region  were  iron-reducing,  sulfate-reducing,  and 
methanogenic.  In  addition,  the  organic  carbon-rich  wetland  supplied  an  abundant 
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quantity  of  electron  donors,  eliminating  a  common  limiting  factor  to  the  reductive 
dechlorination  pathway  (Lorah  and  Olsen,  1999). 

While  Lorah  and  Olsen  (1999)  recognized  how  anaerobic  processes  could 
naturally  biodegrade  chlorinated  hydrocarbons  in  a  wetland,  Amon  et  al.  (2007)  realized 
that  a  combination  of  anaerobic  and  aerobic  processes  in  a  constructed  wetland  could  be 
used  to  biodegrade  chlorinated  ethenes.  Specifically,  aerobic  microbial  populations, 
using  methane  as  a  primary  substrate,  produce  non-specific  oxygenases  which 
cometabolically  degrade  chlorinated  ethenes  like  TCE,  cDCE,  and  vinyl  chloride 
(Bradley,  2003).  An  upward  VFTW,  in  which  the  contaminated  water  flows  first  through 
an  anaerobic  region  and  then  through  an  aerobic  root  zone,  would  allow  both  anaerobic 
and  aerobic  processes  to  operate  with  the  added  benefit  that  methane  produced  in  the 
anaerobic  region  could  be  used  as  the  primary  substrate  to  support  cometabolic 
degradation  of  the  chlorinated  ethenes  in  the  aerobic  region. 

2.2  Frequency  of  Use  of  Upward  VFTW 

Upward  VFTWs  are  a  relatively  recent  development,  and  this  is  reflected  in  the 
paucity  of  information  in  the  literature  relative  to  horizontal  flow  and  downward  VFTWs. 
In  1994,  Hans  Blix  published  an  excellent  survey  of  the  history  of  the  use  of  treatment 
wetlands  and  the  (then)  current  state  of  technology,  covering  104  subsurface  flow  and  70 
FWS  treatment  wetlands  (Blix,  1994).  He  refers  to  two  pilot  scale  uVFTWs  constructed 
in  Australia  for  the  treatment  of  municipal  wastewater.  Both  were  gravel  bed  systems  in 
which  adsorption  and  contact  with  the  root  zone  were  the  primary  mechanisms  of 
contaminant  removal  (Blix,  1994). 
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Farahbakhshazad  and  Morrison  (2000)  profiled  a  uVFTW  built  in  1997  near 
Piracicaba,  Brazil  used  to  remove  ammonia,  nitrate,  and  phosphate  from  the  wastewater 
of  a  small  (1,000  people)  community.  By  obtaining  samples  from  three  locations  within 
the  cell  and  two  depths  per  location,  the  authors  found  that  soil  efficiency  for  the  removal 
of  ammonia  was  not  uniform  throughout  the  cell.  They  speculated  this  was  due  to 
channeling  and  short  circuiting  in  the  subsurface,  but  since  the  authors  lacked  data  on 
hydraulic  head  and  conductivity,  they  were  not  able  to  make  a  definitive  conclusion 
(Farahbakhshazad  and  Morrison,  2000). 

Kassenga  and  colleagues  conducted  bench-scale  experiments  using  an  upflow 
column  designed  to  determine  optimum  soil  composition  and  conditions  for  attenuation 
of  chlorinated  volatile  organic  compounds,  specifically  TCE,  cDCE,  and  1,2-DCA 
(Kassenga  et  al.,  2003).  The  paper,  however,  does  not  refer  to  the  Wright  Patterson 
uVFTW,  which  was  still  very  new  at  the  time  of  the  study.  Kassenga  et  al.  (2003)  state 
“it  was  necessary  to  conduct  bench-scale  studies  to  investigate  the  feasibility  of  using 
constructed  wetlands  for  treatment  of  VOCs,”  implying  that  the  existence  of  a  full-scale 
system  was  not  known.  Pardue  (2005)  published  a  summary  of  bench-scale  research  into 
using  upflow  columns  to  treat  VOC-contaminated  groundwater.  Determinations  were 
made  regarding  the  genetic  identification  of  halorespirers,  physical  and  chemical 
properties  of  various  soil  types,  and  locations  and  rates  of  reactions.  The  effect  on 
treatment  efficiency  was  evaluated  with  respect  to  vegetation  type  and  salinity  and 
inoculation  of  soil  with  halorespirers.  He  refers  to  three  pilot-scale  systems  currently  in 
different  stages  of  EPA  approval:  Superfund  sites  in  Connecticut  and  Massachusetts,  and 


14 


a  site  along  Chesapeake  Bay  within  the  Aberdeen  Proving  Grounds  in  Maryland  (Pardue, 
2005). 

Wallace  and  Kadlec  (2005)  profiled  a  pilot-scale  uVFTW  used  to  degrade  BTEX 
at  the  site  of  a  former  refinery.  Contaminated  water  was  introduced  in  the  bottom  layer, 
then  flowed  vertically  through  gravel  and  sand  layers.  The  water  then  flowed 
horizontally  across  the  sand  layer  to  the  outlet.  This  flow  design  was  noted  as  potentially 
unstable,  and  short-circuiting  was  observed.  In  the  full-scale  system,  the  vertical  upflow 
design  was  abandoned  in  favor  of  a  center-fed,  radial  horizontal  subsurface  flow  design. 

In  a  recently  published  textbook  on  treatment  wetlands  (Kadlec  and  Wallace, 
2009),  the  authors  refer  to  uVFTWs  only  in  passing.  They  make  a  general  statement  that 
uVFTWs  are  used  to  promote  anaerobic  degradation.  In  the  chapter  on  organic  chemicals 
the  uVFTW  at  Wright-Patterson  AFB  is  specifically  referred  to  in  only  two  sentences  and 
no  information  is  given  as  to  how  the  treatment  wetland  operates.  The  authors  cite  two 
AFIT  Master’s  theses  in  this  passage,  although  interestingly,  the  two  cited  theses 
(Entingh,  2002;  and  Blalock,  2003)  do  not  deal  with  degradation  of  organic  chemicals, 
but  rather  with  characterizing  the  hydraulic  properties  of  the  treatment  wetland  cell.  In 
any  event,  with  the  exception  of  the  Wright-Patterson  AFB  uVFTW,  Kadlec  and  Wallace 
(2009)  do  not  cite  any  examples  of  the  use  of  uVFTWs  to  treat  chlorinated  hydrocarbons. 

The  first  peer-reviewed  article  on  the  Wright  Patterson  AFB  uVFTW  was 
published  by  Amon  et  al.  (2007)  of  Wright  State  University.  The  authors  were  aware  of 
no  other  uVFTWs  then  in  existence  being  used  to  treat  chlorinated  organic  compounds. 
See  section  2.7  Case  Studies  for  further  information  on  this  paper,  as  well  as  master’s  and 
PhD  research  involving  the  uVFTW  at  WPAFB. 
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2.3  Typical  Values  for  Environmental  and  Constructed  Wetland  Design 

Parameters 

Lorah  and  Olsen  (1999)  determined  that  the  wetland  sediment  thought  responsible 
for  the  reductive  dechlorination  of  VOCs  was  1.8-3. 6m  thick  and  comprised  of  two 
layers;  a  lower  silty  clay  layer  low  in  organic  carbon  (1%),  and  an  upper  peat  and  clay 
layer  high  in  organic  carbon  (18%).  Average  linear  groundwater  flow  velocity  in  the 
wetland  sediment  was  estimated  at  0.6  m  yr'1.  The  redox  conditions  in  the  wetland 
sediment  were  extensively  characterized,  and  the  authors  noted  that  the  reductive 
dechlorination  pathway  was  most  effective  under  methanogenic  conditions.  The  authors 
did  not  evaluate  the  porosity  or  hydraulic  conductivity  of  the  wetland  sediments  (Lorah 
and  Olsen,  1999). 

In  the  study  of  Kassenga  et  al.  (2003)  discussed  earlier,  two  soil  mixtures  were 
compared:  one  containing  20%  sand  and  80%  peat,  and  the  other  containing  20%  sand, 
40%  peat,  and  40%  Bion  soil,  a  commercially  available  product  derived  from  agricultural 
waste  (Kassenga  et  al.,  2003).  A  summary  of  soil  parameters  is  presented  in  Table  1. 
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Table  1:  Soil! 

Parameters  from  Kassenga  et  al.  (2003) 

Parameter 

20%  sand 

80%  peat 

20%  sand 

40%  peat  and  40%  Bion  soil 

Porosity  (n) 

0.76+0.04 

0.72+0.06 

Retardation  Factor  (R) 

cDCE  3.51 

TCE  2.75 

1,1,1-TCA  3.69 

VC  1.07 

cDCE  7.32 

TCE  4.64 

1,1,1-TCA  6.56 

VC  1.20 

Hydraulic  conductivity  (K) 

0.43  m  d1 

0.11  md'1 

Critical  Flow 

Total  flow  for  18m  x  36m-cell 

0.080  m3  d'1  m"2 

51.8  m3  d'1 

0.074  m3  d"1  m  “ 

48.0  m3  d'1 

Minimum  bed  depth  for 
effective  VC  removal 

9.00  m 

2.10m 

Pardue  (2005)  did  microcosm  and  mesocosm  work  similar  in  nature  to  Kassenga 
et  al.  (2003)  (Kassenga  did  PhD  work  under  Pardue).  Pardue  (2005)  calculated  removal 
rates  for  1,1, 2, 2  tetrachloroacetate,  TCE,  cDCE,  and  1,2  dichloroacetate  and  determined 
that  most  removal  occurred  in  the  bottom  10  cm  of  the  column  (Pardue,  2005).  This  was 
based  on  a  flow  of  6.41E-03  gallons  per  minute  per  nT.  Applying  this  value  to  the  648 
m“  treatment  wetland  at  WPAFB  give  4.2  gallons  per  minute,  which  is  slightly  less  than 
the  5.5  gallons  per  minute  reported  by  Corbin  (2008). 

The  soil  type  used  in  the  Wright  Patterson  uVFTW  was  characterized  as  Westland 
soil  with  silt  and  clay  inclusions,  a  soil  type  common  to  local  area  fens  (Amon  et  al., 
2007).  A  fen  is  a  peat-forming  wetland  fed  by  groundwater  rather  than  precipitation, 
which  results  in  a  generally  high  nutrient  level  in  the  soil  (U.S.  EPA,  2009).  Because  the 
direction  of  flow  is  toward  the  surface  in  a  fen,  this  type  of  soil  should  be  well-suited  to 
use  in  uVFTW  (Table  2). 
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Table  2:  Selected  Soil  Data  from  Amon  et  al.  (2007) 


Effective  Porosity 

26% 

Hydraulic  conductivity 

0.81  md'1 

Inflow  rate 

54.7  mJ  d'1 

Depth  of  treatment  wetland 

1.45  m 

2.4  Modeling  Biodegradation  Kinetics  of  Chlorinated  Ethenes  in  Constructed 

Wetlands 

The  choice  of  kinetics  used  to  model  biodegradation  of  chlorinated  ethenes  in  a 
treatment  wetland  cannot  be  made  without  an  understanding  of  the  primary  degradation 
processes  occurring  in  the  wetland  (Bradley,  2003).  The  primary  processes  that  occur  in 
a  uVFTW  are: 

1)  Anaerobic  reductive  dechlorination 

2)  Aerobic  cometabolism  (methane  as  primary  substrate) 

In  anaerobic  reductive  dechlorination,  microorganisms  use  Hi  from  fermented 
organic  compounds  present  in  organic  rich  soil  as  the  electron  donor,  and  PCE  or  TCE  as 
the  terminal  electron  acceptor  (Figure  4).  This  process,  also  known  as  chlororespiration, 
results  in  a  net  energy  gain  for  the  microorganisms.  Anaerobic  reductive  dechlorination 
results  in  a  chlorine  ion  being  stripped  off  the  PCE  molecule  and  replaced  with  a 
hydrogen  ion. 
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Figure  4:  Anaerobic  Reductive  Dechlorination  of  PCE  to  TCE  and  TCE 

to  cDCE 
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Electron  acceptor 

-PCE,  TCE  primarily  (DCE  and  VC  to  much  less  extent) 

-Methanogenic  conditions  provide  best  rates  of  reductive  dechlorination. 


Source:  material  adapted  from  Vogel  et  al.  (1987) 


Anaerobic  reductive  dechlorination  is  a  stepwise  process,  meaning  PCE  is  first 
degraded  to  TCE,  which  is  then  degraded  to  DCE.  Here,  the  process  tends  to  stall  due  to 
the  inhibitory  effect  of  increasing  concentrations  of  daughter  products,  as  well  as  the 
increasingly  reducing  conditions  required  to  strip  additional  chlorine  ions.  As  chlorinated 
ethenes  become  reduced,  the  effectiveness  of  reducing  pathways  is  less. 

Microbes  capable  of  using  chlorinated  ethenes  as  electron  acceptor  must  compete 
for  electron  donors  with  other  microorganisms  and  other  redox  processes.  The  redox 
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potential  governs  how  successful  the  process  of  anaerobic  reductive  dechlorination  will 


be  (Table  3). 


Table  3:  Gibbs  Free  Energy  Change  for  Sel 

lected  Processes 

Redox  Environment 

Gibbs  free  energy  change 

Oxygen  reducing 

-502  kJ/mol 

Nitrate  reducing 

-477  kJ/mol 

Manganese  reducing 

-339  kJ/mol 

-164  kJ/mol 

PCE  to  TCE1 

-161  kJ/mol 

TCE  to  cDCE1 

-141  kJ/mol 

cDCE  to  VC1 

-154  kJ/mol 

VC  to  Ethene1 

Iron  reducing 

-117  kJ/mol 

Sulfate  reducing 

-105  kJ/mol 

CO2  reducing 
(methanogenic) 

-92  kJ/mol 

Source:  (U.S.  EPA,  1999), 
except  for:  1  (He  et  al.,  2002) 


Table  3  indicates  that  anaerobic  reductive  dechlorination  will  occur  to  a  much 
reduced  extent  under  nitrate  and  manganese  reducing  conditions  as  compared  to 
methanogenic  conditions.  Bradley  (2003)  noted  that  PCE  will  undergo  anaerobic 
reductive  dechlorination  under  all  but  aerobic  conditions,  TCE  will  undergo  anaerobic 
reductive  dechlorination  under  iron  reducing,  sulfate  reducing,  or  methanogenic 
conditions,  DCE  will  undergo  anaerobic  reductive  dechlorination  under  either  sulfate 
reducing  or  methanogenic  conditions,  and  VC  will  undergo  anaerobic  reductive 
dechlorination  only  under  methanogenic  conditions  (Bradley,  2003). 
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A  primary  substrate  to  support  aerobic  cometabolism  of  chlorinated  ethenes  (i.e., 
methane)  may  be  produced  from  COo  and  Hi  via  anaerobic  methano genesis.  This 
methane  is  then  used  by  methylotrophs,  which  produce  methane  monooxygenase  (MMO) 
enzyme  under  aerobic  conditions  to  catalyze  oxidation  of  methane  to  produce  CCK 
MMO  is  a  nonspecific  enzyme  that  is  also  capable  of  aerobically  degrading  TCE,  DCE, 
and  VC  to  CO2,  without  accumulations  of  harmful  intermediate  products.  The 
methylotroph  itself  oxidizes  the  primary  substrate  for  growth  and  energy,  but  appears  to 
derive  no  benefit  from  the  cometabolic  oxidation  of  the  chlorinated  ethene  (Figure  5). 

The  extent  of  cometabolism  is  typically  limited  by  the  availability  of  the  primary 
substrate,  methane  (Bradley,  2003). 


Figure  5:  Aerobic  Oxidation  of  TCE  to  CO2  (cometabolism) 


Source:  Material  adapted  from  McCarty  (1997)  in  Subsurface  Restoration,  p374 
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With  an  understanding  of  the  primary  degradation  processes  that  occur  in  a 
uVFTW,  one  can  then  look  at  kinetics.  Monod  kinetics  relates  the  growth  rate  of  a 
microorganism  to  the  concentration  of  a  limiting  substrate,  and  it  has  been  used  to 
describe  the  growth  rate  of  microorganisms  engaged  in  anaerobic  reductive 
dechlorination  (Suarez  and  Rifai,  1999). 


5 

Lhnax 


(10) 


where, 


p:  growth  rate  of  biomass  (microorganism)  [T1] 

Pmax-  max  growth  rate  of  biomass  (microorganism)  [T  1  ] 

S:  growth-limiting  substrate  concentration  [M  L  3] 

Ks:  substrate  concentration  that  results  in  a  growth  rate,  /4  of  p  max/2 

[M  L'3] 


The  rate  of  biomass  growth  is  related  to  the  rate  of  substrate  consumption  by 
equation  (11)  (Kovarova-Kovar  and  Egli,  1998). 


Y  dS 
u  =  —  x  — 
B  dt 


(11) 


where, 

Y:  yield  coefficient,  mass  biomass  per  mass  substrate  [M  M"1] 

B:  biomass  concentration  [M  IT  ] 
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Substituting  (11)  into  (10),  rearranging  the  expression,  and  introducing  a  negative 
sign  to  reflect  the  fact  that  substrate  is  being  consumed  gives  an  expression  for  the  rate  of 
change  in  substrate  (chlorinated  ethene)  concentration. 


dS 

dt 


(12) 


where, 

Y:  yield  coefficient,  mass  biomass  per  mass  substrate  [M  M ’] 

B:  biomass  concentration  [M  L‘ ] 

From  Equation  (12),  it  can  be  seen  that  when  S«Ks: 


Rearranging, 


dS 

dt 


Bmax 

Y 


B 


(13) 


dS  —  /  Bmax  „ 

~dt  ~  V  Y  x  Ks  ) 


It  can  be  seen  from  this  that  when  S«Ks,  dS/dt  can  be  approximated  by  the  first  order 
expression: 


dS 

dt 


=  -kS 


(15) 


where, 
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(16) 


ftmax  ^ 

.  YxKs  ) 


First-order  rate  constants  that  have  appeared  in  the  literature  are  reported  in  Table 
4  and  Table  5. 


Table  4:  First  Order  Degradation  Rate  Constants 


PCE  to  TCE 

TCE  to  cDCE 

cDCE  to  VC 

VC  to  Ethene 

Anaerobic 

reductive 

dechlorination 

Mean:  .010 
90th%:  .022 
Range:  0-0.080 
(Suarez  and 
Rifai,  1999) 

Mean:  .003 
90th%:  .005 
Range:  0-0.023 
(Suarez  and 

Rifai,  1999) 

Mean:  .003 
90th%:  .005 
Range:  .001- 
.006 

(Suarez  and 
Rifai,  1999) 

Mean:  .001 
90th%:  n/a 

Range:  0-0.007 
(Suarez  and 

Rifai,  1999) 

All  rates  listed  are  in  units  of  [d  ].  Data  were  taken  from  Table  7. 


Table  5:  First  Order  Degradation  Rate  Constants,  Aerobic  Cometabolism 


PCE 

TCE 

cDCE 

VC 

Aerobic  \  / 

Mean:  .15+.02 

Mean:  .59+.07 

cometabolism  \  / 

(Powell  et  al., 

(Powell  et  al., 

\  / 

2010) 

2010) 

Mean:  .948 

Mean:  .720 

Mean:  1.730 

y 

90th%:  - 

90th%:  1.012 

90th%:  - 

/  \ 

Range:  .105- 

Range:  .390- 

Range:  1.50- 

/  \ 

1.41 

1.15 

1.96 

/  \ 

(Suarez  and 

(Suarez  and 

(Suarez  and 

/  \ 

Rifai,  1999) 

Rifai,  1999) 

Rifai,  1999) 

All  rates  listed  are  in  units  of  [d  ].  Data  from  Suarez  and  Rifai  (1999)  were  taken  from 
Table  7.  Rates  from  Powell  et  al.  (2010)  are  pseudo  first-order,  and  are  from  lab  studies. 
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2.5  Relationship  Between  Hydraulic  Conductivity,  Residence  Time,  and 

Treatment  Efficiency 

The  significance  of  the  problem  of  hydraulic  short-circuiting  is  difficult  to 
determine  in  the  case  of  uVFTWs  due  to  the  fact  that  uVFTWs  are  relatively  recent 
innovations  and  not  many  are  in  use.  Where  they  have  been  used,  hydraulic  short- 
circuiting  is  noted  as  a  potential  problem,  but  the  degree  to  which  it  is  a  problem  is  not 
quantified  (Farahbakhshazad  and  Morrison,  2000;  Wallace  and  Kadlec,  2005). 

Hydraulic  short-circuiting  has  been  studied  much  more  in  HSFW  treatment 
wetlands  and  dVFTWs.  Araujo  et  al.  (2008)  modeled  the  problem  of  clogging  in  HSFW 
treatment  wetlands,  noting  that  dead  zones  and  hydraulic  short-circuiting  are  common 
ramifications  of  clogging  (Araujo  et  al.,  2008). 

Chazarenc  et  al.  (2003)  noted  hydraulic  short-circuiting  and  dead  zones  in 
conducting  pulse  tracer  experiments  to  determine  hydraulic  residence  time  distributions 
in  a  HSFW  in  France.  He  attributed  this  to  the  structural  design  of  the  bed,  where 
influent  is  injected  at  one  point  in  the  comer  of  the  bed  and  effluent  is  removed  at  the 
opposite  comer,  a  design  which  promotes  hydraulic  short-circuiting  (Chazarenc  et  al., 
2003). 

Dittrich  (2006),  in  a  review  of  Hungarian  usage  of  HSFWs  and  dVFTWs,  noted 
that  hydraulic  short-circuiting  was  a  significant  problem  in  both  types  due  in  part  to 
distribution  pipe  construction  deficiencies  and  the  shape  of  the  wetland  (Dittrich,  2006). 
He  proposed  the  use  of  only  gravel  as  filter  media  in  HSFWs  in  order  to  avoid  hydraulic 
short-circuiting,  and  avoiding  the  use  of  gravel  as  filter  media  in  dVFTWs  in  order  to 
increase  hydraulic  residence  time  (Dittrich,  2006). 
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Lightbody  et  al.  (2008)  observed  hydraulic  short-circuiting  in  a  three-cell  FWS 
wetland  in  Georgia.  Between  three  and  six  fast  flow  paths  were  observed  which  carried 
water  at  least  10  times  as  fast  as  other  areas  of  the  cell,  and  20-70%  of  the  flow  had  a 
residence  time  less  than  the  nominal  residence  time  (Lightbody  et  al.,  2008).  The  authors 
also  note  that  the  wetland  was  carefully  designed  to  ensure  uniform  flow,  with  careful 
planting  of  vegetation,  multiple  inlet  and  outlet  structures,  and  a  flat  base  to  the  cell 
(Lightbody  et  al.,  2008).  The  fact  that  short-circuiting  still  occurred  may  indicate  the 
difficulty  in  eliminating  this  problem  by  design. 

In  the  examples  cited  above,  hydraulic  short-circuiting  is  determined  by 
calculating  the  nominal  or  theoretical  residence  time  by  dividing  water  volume  in  the  cell 
by  the  volumetric  flow  rate.  This  is  then  compared  to  the  mean  or  first  moment  of  the 
residence  time  distribution  function  (RTDF).  If  the  first  moment  of  the  RTDF  is  less  than 
tbar,  hydraulic  short-circuiting  is  indicated.  Figure  6  illustrates  an  example  where  there 
appears  to  be  short-circuiting,  and  flow  is  not  going  through  the  entire  cell  volume. 
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Figure  6:  Hypothetical  RTDF 


0.45 


Typical  treatment  efficiencies  vary  widely,  due  to  wetland  design,  type  of 
contaminant,  and  hydraulic  retention  time.  Farahbakhshazad  and  Morrison  (2000) 
measured  removal  rates  for  nitrate,  phosphate,  and  ammonia  for  a  uVFTW  in  Brazil. 
Concentration  removal  for  phosphate  was  93%,  78%  for  nitrate,  and  50%  for  ammonia 
(Farahbakhshazad  and  Morrison,  2000).  The  total  depth  of  the  bed  was  0.45  m,  and  the 
upper  0.25  m  was  vermiculite  with  about  35-60%  clay  content.  Sampling  was  conducted 
at  three  sites  within  the  uVFTW,  at  two  depths  per  site  -  0.25  m  (bottom  of  the 
vermiculite  layer)  and  0.35  m  (bottom  gravel  layer).  The  mechanism  of  removal  was 
adsorption  in  the  bed  and  precipitation  at  the  surface  for  phosphate,  and  plant  uptake  for 
nitrate  and  ammonia  (Farahbakhshazad  and  Morrison,  2000). 
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Kassenga  et  al.  (2003),  in  column  experiments  of  cDCE  and  VC  degradation, 
determined  the  minimum  bed  depth  to  degrade  influent  concentrations  of  both 
contaminants  of  13.67  mg/L  to  below  the  maximum  contaminant  levels  (MCLs)  specified 
in  National  Primary  Drinking  Water  Regulations.  VC  was  determined  to  be  the  limiting 
contaminant,  requiring  the  deepest  bed  depth  to  degrade  the  contaminant  to  below  its 
MCL  of  2  pg/L.  Based  on  a  first-order  rate  constant  of  0.56  day"1,  a  bed  depth  of  2.1m 
was  required  (Kassenga  et  al.,  2003). 

Lee  et  al.  (2004)  evaluated  the  removal  efficiency  of  a  HSFW  designed  to  treat 
swine  waste  effluent.  Removal  efficiency  for  a  variety  of  effluent  constituents  was 
evaluated  for  three  hydraulic  residence  times:  4.3  days,  8.5  days,  and  14.7  days  (Figure 
7).  There  was  no  consistent  trend  in  removal  efficiency  as  a  function  of  hydraulic 
residence  time.  Several  constituents  (NO3-N,  total  phosphorus,  and  PO4-P)  exhibited  a 
positive  correlation  between  removal  efficiency  and  hydraulic  residence  time,  as  would 
be  expected,  while  others  (COD,  NO2-N,  NH4-N,  and  total  nitrogen)  appeared  to  have  an 
optimum  removal  efficiency  when  the  hydraulic  residence  time  was  8.5  days.  The 
authors  note  that  the  profiles  for  nitrate,  nitrite,  and  ammonia  may  be  due  to  the  poor 
nitrification  capability  of  the  treatment  wetland  (Lee  et  al.,  2004). 
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Figure  7;  Percent  Removal  as  a  Function  of  Hydraulic  Residence  Time 


- COD 

- BOD5 

- NH-N 

- N02-N 

- N03-N 

TKN 

- TN 

P04-P 


Source:  Adapted  from  Table  1  and  Table  2  of  Lee  et  al.  (2004) 


Dierberg  et  al.  (2005)  conducted  a  rhodamine  dye  tracer  study  on  a  FWS  in  south 
Florida  that  exhibited  significant  hydraulic  short-circuiting.  They  determined  that  a  small 
volume  of  the  wetland  conveyed  44%  of  the  volumetric  flow,  and  that  this  reduced  the 
treatment  efficiency  for  phosphorus  in  two  ways.  First,  through  the  shorter  residence 
time;  and  second,  through  impaired  physical  processes  such  as  sedimentation  and  plant 
uptake.  First-order  degradation  rate  constants  for  removal  of  phosphorus  in  the  high 
residence  time  zones  were  twice  those  of  the  short-circuit  zones  (0.50  d_1  vs.  0.24  d ') 
(Dierberg  et  al.,  2005). 
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As  noted  above,  although  hydraulic  short-circuiting  of  constructed  wetlands 
appears  to  be  a  problem,  proposed  engineering  solutions  to  the  problem  and  attempts  to 
simulate  the  effects  of  implementing  proposed  solutions  on  wetland  performance  are  only 
found  in  the  literature  infrequently.  In  one  example,  Chazarenc  et  al.  (2003) 
recommended  “wide  centered  injection”  at  the  influent  of  HSFWs  in  order  to  promote 
mixing  and  prevent  dead  zones,  but  did  not  offer  a  more  detailed  explanation,  and  made 
no  attempt  to  model  such  a  design.  Lightbody  et  al.  (2009)  modeled  a  short-circuiting 
FWS  treatment  wetland  using  a  stream  tube  model  with  dispersion  and  incorporating 
transverse  deep  zones  as  an  engineering  solution  intended  to  mitigate  short-circuiting. 

The  deep  zones  stretched  the  entire  width  of  the  treatment  wetland,  transverse  to  the 
direction  of  flow,  and  were  vegetation-free  pools  from  bottom  to  top.  The  solution  was 
thought  to  improve  treatment  efficiency  via  two  mechanisms.  First,  the  water  is  rapidly 
mixed  in  a  lateral  direction,  which  dilutes  contaminant  concentrations.  Second, 
transverse  deep  zones  disrupt  fast  flow  paths  and  reduce  the  likelihood  of  any  one  path 
extending  the  entire  length  of  the  treatment  cell  (Lightbody  et  ah,  2009).  The  authors 
founds  that  transverse  deep  zones  would  increase  contaminant  removal,  provided  the 
zone  is  long  enough  (in  the  direction  of  flow)  to  completely  dissipate  the  momentum  of 
incoming  fast-path  water.  Further,  contaminant  removal  increases  with  the  number  of 
deep  zones  up  to  a  maximum,  beyond  which  wetland  area  reserved  for  contaminant 
treatment  drops  too  low  and  contaminant  removal  goes  down  (Lightbody  et  al.,  2009) 
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2.6  Use  of  Models  to  Predict  Treatment  Efficiency  of  Constructed  Wetlands 

Lightbody  et  al.  (2009)  showed  how  modeling  could  be  useful  in  predicting 
efficiency  and  assisting  in  wetland  design.  This  section  will  examine  other  such  uses  of 
modeling.  As  the  understanding  of  the  hydrological,  physical,  and  chemical  processes 
that  occur  in  treatment  wetlands  has  improved,  so  has  the  capability  to  model  these 
systems.  Models  of  treatment  wetlands  were  originally  conceived  as  “black  boxes,” 
where  concentration  in  was  linked  to  concentration  out  by  first-order  decay  coefficients 
(Langergraber,  2008).  Most  models  in  use  today  are  deterministic  mathematical  models 
which  use  numerical  methods  to  arrive  at  solutions  to  the  equations  of  flow  and  transport. 
These  numerical  models  allow  for  heterogeneity  in  space  (e.g.,  spatially  varying 
hydraulic  conductivity)  and  variations  in  time  (e.g.,  temporally  varying  influent 
concentrations)  (Konikow  et  al.,  2007). 

Langergraber  (2008)  provides  a  review  of  mechanistic  models  divided  into  three 
categories:  models  of  hydraulic  behavior  and  transport  only;  models  of  reactive  transport 
under  saturated  conditions;  and  models  of  reactive  transport  under  variably  saturated 
conditions.  Of  these,  the  models  of  reactive  transport  in  saturated  conditions  are  relevant 
to  this  research.  Using  a  mixing  cell  method,  a  HFTW  is  divided  into  numerous  equal 
sized  cells,  within  which  degradation  occurs  by  first  order  decay.  This  method  is  a 
simplification  of  the  advection-dispersion  equation,  and  provides  a  better  fit  than  plug 
flow.  A  mechanistic  model  using  STELLA  was  also  presented.  The  model  incorporates 
six  state  variables  -  carbon  cycle,  nitrogen  cycle,  oxygen  balance,  bacteria  growth 
(heterotrophic  and  autotrophic),  and  water  budget.  It  predicted  effluent  BOD,  organic 
nitrogen,  ammonium,  and  nitrate  concentrations  well.  STELLA  was  also  used  in  another 
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approach  to  model  nitrogen  transformation  in  HFTWs.  Other  models  use  a  variety  of 
combinations  of  continuously  stirred  tank  reactors  (CSTRs)  and  plug  flow  reactors  with 
the  goal  of  simplifying  the  model  by  reducing  the  number  of  parameters  required 
(Langergraber,  2008). 

Langergraber  et  al.  (2009)  demonstrates  the  multi-component  reactive  transport 
module  CW2D,  designed  to  model  the  degradation  of  municipal  wastewater  in  a  HFTW. 
CW2D  is  capable  of  modeling  12  components  and  9  processes  using  Monod  kinetics  to 
describe  degradation.  CW2D  is  imbedded  within  the  HYDRUS  transport  program  for 
variably  unsaturated  flow  (Langergraber  et  al.,  2009). 

One  of  the  most  widely  used  groundwater  modeling  software  packages  is 
MODFLOW,  developed  by  the  U.S.  Geological  Survey  (Konikow  et  al.,  2007). 
MODFLOW  uses  the  finite  difference  method  to  arrive  at  numerical  solutions  to  the 
groundwater  flow  equation  (Harbaugh,  2005).  MODFLOW  is  a  modular  software 
package  which  can  be  coupled  with  other  packages,  particularly  packages  that  model 
transport.  RT3D,  Reactive  Transport  in  3  Dimensions,  is  a  separate  software  package 
that  simulates  transport  and  degradation  (Clement,  1997).  RT3D  uses  the  output  of 
MODFLOW  to  determine  advective  velocities  to  input  into  the  transport  equations. 

RT3D  contains  built-in  programs  to  simulate  many  different  fate  and  transport  processes, 
including  a  model  for  both  anaerobic  and  aerobic  degradation  of  PCE  and  TCE  (Clement, 
1997).  Although  MODFLOW  and  RT3D  can  be  used  directly,  it  is  also  common  to  use  a 
more  user-friendly  interface  program  such  as  the  Groundwater  Modeling  System  (GMS). 
GMS  is  a  Windows-based  interface  to  MODFLOW-2005  and  RT3D  created  by  the 
Department  of  Defense  that  integrates  multiple  programs  for  subsurface  flow  and 
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contaminant  fate  and  transport  (U.S.  A.C.E.,  2010).  This  research  effort  will  use 
MODFLOW-2005  and  RT3D  within  the  GMS  interface,  so  further  information  about 
them  is  included  in  chapter  3. 

2.7  Case  study:  uVFTW  at  WPAFB 

This  section  will  be  devoted  to  summarizing  prior  research  into  the  hydraulic 
behavior  and  degradation  characteristics  of  the  uVFTW  at  WPAFB.  The  first  peer- 
reviewed  article  on  the  Wright  Patterson  AFB  uVFTW  was  published  in  2007  (Amon  et 
al.,  2007).  Amon  et  al.  (2007)  gives  a  detailed  account  of  construction  and  operation  of 
the  uVFTW,  to  include  inflow  rate,  physical  dimensions  of  the  cell,  composition  and 
depths  of  the  different  layers  of  the  cell,  influent  piping,  types  of  vegetation  planted,  and 
piezometer  layout  (Amon  et  al.,  2007).  Amon  notes  that  while  anaerobic  reductive 
dechlorination  of  PCE  is  a  major  process  in  the  wetland,  at  shallow  depths,  due  to  the 
availability  of  methane  and  root-transported  oxygen,  oxidative  cometabolism  is  an 
important  degradation  process,  as  well  (Amon  et  al.,  2007).  Amon  et  al.  (2007)  suggest 
the  inclusion  of  wood  chips  throughout  the  cell,  as  well  as  filling  the  cell  with  dirt  as  it 
fills  with  water  (during  construction),  as  means  of  avoiding  uneven  compaction  and 
hydraulic  short-circuiting  (Amon  et  al.,  2007). 

Many  master’s  and  PhD  students  at  Wright  State  University  and  the  Air  Force 
Institute  of  Technology  have  done  research  on  the  uVFTW  at  WPAFB.  Mesocosm 
studies  of  the  soil  type  used  in  the  uVFTW  include  Tritschler  (2007),  Gruner  (2008),  and 
Powell  (201 1).  Studies  examining  biodegradation  characterization  of  the  treatment 
wetland  include  Opperman  (2002),  Clemmer  (2003),  Kovacic  (2003),  Sobolewski  (2004), 
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Waldron  (2007),  and  Thompson  (2008).  Research  into  the  hydraulics  of  the  cell  includes 
Entingh  (2002),  Blalock  (2003),  and  Corbin  (2008). 

Tritschler  et  al.  (2007)  conducted  mesocosm  experiments  to  investigate  the  effects 
of  vegetation  and  seasonal  changes  on  redox  parameters  within  the  mesocosm.  She 
determined  that  selection  of  plant  species  had  a  large  impact  on  the  redox  conditions 
(Tritschler  et  al.,  2007).  Therefore,  the  ability  of  a  treatment  wetland  to  degrade 
contaminants  would  be  sensitive  to  selection  of  plant  species  (Tritschler  et  al.,  2007). 
Gruner  (2008)  conducted  mesocosm  experiments  to  characterize  bacterial  populations  of 
methanogens,  methanotrophs,  and  ammonia  oxidizers  with  respect  to  planted  versus  non- 
planted  mesocosms.  Powell  (201 1)  conducted  field  work  characterizing  contaminant 
concentrations  in  the  shallow  vegetated  zone,  as  well  as  laboratory  microcosm  work  to 
determine  degradation  rates  for  TCE,  cDCE,  and  1,1,1-TCA  in  the  presence  of  methane 
oxidizers.  TCE  and  cDCE  were  found  to  degrade  with  first-order  degradation  rates, 
while  1,1,1-TCA  was  not  found  to  degrade  (Powell  et  al.,  2010). 

Kovacic  (2003)  conducted  field  sampling  to  characterize  organic  acid  and 
inorganic  anion  concentrations  in  the  wetland  cell.  He  found  that  organic  acid 
concentrations  decreased  by  93%  over  the  11  months  between  December  2002  and 
January  2003  (Kovacic,  2003).  From  measurements  of  dissolved  oxygen  and  oxidative- 
reductive  potential,  he  determined  that  the  base  of  the  wetland  was  aerobic,  but  that 
quickly  changed  to  anaerobic  in  the  layer  immediately  above  the  base.  He  also  found  that 
nitrate  and  sulfate  reducing  conditions  were  correlated  with  higher  concentrations  of 
lactate  and  formate  (Kovacic,  2003). 
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Waldron  (2007)  conducted  field  sampling  research  to  characterize  concentration 
profiles  for  PCE,  TCE,  trans-DCE,  cDCE,  VC,  and  ethane  at  each  of  the  three  depths  of 
the  wetland  (Waldron,  2007).  This  work  was  a  continuation  of  similar  field  work  done 
by  Sobolewski  (2004),  Clemmer  (2003),  and  Opperman  (2002).  Waldron  presented 
summary  concentration  data  for  all  the  sampling  -  December  2001,  January  2003,  Fall 
2003,  September  2006,  October  2006,  November  2006,  and  December  2006  (Waldron, 
2007).  This  allowed  him  to  state  some  observations  regarding  changes  in  the  uVFTW’s 
treatment  efficiency  over  time: 

1)  Influent  PCE  concentration  increased  by  94%  from  Fall  2003  to  Fall  2006,  but 
treatment  efficiency  (measured  at  the  top  layer,  not  the  effluent  weir)  held 
constant  at  99%.  This  has  implications  for  the  maximum  influent  contaminant 
level  the  wetland  cell  can  handle  without  the  effluent  exceeding  maximum 
contaminant  levels  (MCLs). 

2)  Reduction  of  PCE  concentrations  was  less  measured  at  the  effluent  weir  (85.4%) 
compared  to  the  top  layer  of  the  wetland  cell  (99%).  Waldron  concluded  from 
this  that  some  influent  water  was  bypassing  treatment  and  short-circuiting  directly 
to  the  effluent  weir. 

3)  TCE  concentrations  in  the  middle  and  upper  layers  declined  substantially  from 
Fall  2003  to  Fall  2006.  Waldron  attributes  this  to  continued  maturation  of  the 
wetland  cell. 

4)  The  concentration  of  trans-DCE  increased  significantly  over  the  years,  becoming 
the  dominant  DCE  isomer  in  the  wetland  cell. 

5)  The  concentration  of  VC  decreased  from  2003  to  2006  by  approximately  an  order 
of  magnitude,  further  indicating  maturation  of  the  cell  and  somewhat  allaying 
concerns  about  build  up  of  this  end  product. 

Data  presented  by  Waldron  show  classic  reductive  dechlorination  (Figure  8).  As 
PCE  is  degraded  from  the  lower  to  middle  layer,  concentrations  of  daughter  products 
increase.  These  concentrations  drop  off  from  middle  to  upper  layers,  suggesting  either  a 
continuation  of  sequential  anaerobic  reductive  dechlorination,  or  that  aerobic  processes 
become  more  prominent  in  the  upper  layer  (Waldron,  2007).  Not  shown  in  Figure  8  due 
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to  scale  are  concentration  profiles  for  VC  and  Ethane,  which  show  zero  concentration  at 
influent,  rising  to  a  peak  at  the  middle  layer,  then  falling  back  to  almost  zero  at  effluent. 


Figure  8:  Average  Contaminant  Concentration  by  Layer 
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Source:  Waldron  (2007) 


Thompson  (2008)  used  STELLA  to  model  the  movement  of  oxygen  from  roots 
into  the  rhizosphere,  as  well  as  aerobic  cometabolism  of  TCE  (Thompson,  2008).  He 
determined  that  low  organic  carbon  concentration  in  the  influent,  low  copper,  high 
oxygen,  and  high  methane  concentrations  improved  destruction  of  TCE.  Treatment 
efficiency  was  most  affected  by  hydraulic  loading,  dropping  from  80%  to  20%  with  an 
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2-1  2-1 

increase  in  loading  from  1  L  m"“  hr'  to  4  L  m‘~  hr'  (Thompson,  2008).  The  current 
loading  rate  for  the  uVFTW  is  about  30  m3  d'1,  which  is  0.5  L  m  2  hr  1  (Corbin,  2008). 

Entingh  (2002)  was  the  first  to  look  at  the  hydraulic  behavior  of  the  uVFTW.  He 
describes  installation  of  the  grid  of  66  piezometer  nests,  each  of  which  contain  3 
piezometers  screened  in  the  lower,  middle,  and  upper  layers  of  the  wetland  cell  (Entingh, 
2002).  Head  was  measured  at  all  locations,  and  slug  tests  were  conducted  using  the 
Bouwer  and  Rice  (1976)  method  at  all  locations  in  order  to  calculate  values  of  hydraulic 
conductivity.  Contour  plots  of  head  and  hydraulic  conductivity  were  constructed,  which 
indicated  regions  of  preferential  flow.  On-site  evaluation  confirmed  this,  as  areas  of 
higher  head  were  soft  and  could  not  support  the  weight  of  foot  traffic.  One  area  exhibited 
fluidization,  with  significant  mounding  and  water  flowing  out  of  the  ground. 

Because  some  piezometers  were  not  recovering  quickly  enough  to  allow  sample 
collection  to  test  for  PCE,  these  piezometers  were  developed  by  pumping  water  into 
them.  26  piezometers  in  the  top  layer  and  10  in  the  middle  layer  were  initially  developed 
using  this  technique.  This  seems  to  have  altered  soil  properties  immediately  surrounding 
the  developed  piezometers,  resulting  in  changes  in  calculated  hydraulic  conductivity 
values  of  1  or  2  orders  of  magnitude  in  some  cases  (Table  6).  In  order  to  ensure 
conditions  over  the  entire  wetland  were  the  same,  the  remaining  piezometers  in  the  top 
and  middle  layer  were  subsequently  developed  (Entingh,  2002). 
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Table  6:  Geometric  Means  of  Hydraulic  Conductivity  (m  d'1) 


Top  Soil  Layer 

Middle  Soil  Layer 

Bottom  Soil  Layer 

Before  Developing 

0.0011 

0.0066 

14.8 

After  Developing 

0.09 

0.06 

14.8 

Source:  Entingh  (2002) 


Calculated  hydraulic  residence  times  provided  further  evidence  of  hydraulic 
short-circuiting.  These  ranged  from  16.5  hours  to  15  days  (mean  3  days). 

Approximately  64%  of  the  particles  from  a  MODPATH  simulation  had  residence  times 
less  than  the  mean  (Entingh,  2002). 

The  next  year,  Blalock  (2003)  also  examined  the  hydraulic  behavior  of  the 
uVFTW  at  WPAFB.  Whereas  Entingh  (2002)  used  the  Bouwer  and  Rice  (1976)  method 
to  determine  hydraulic  conductivity,  Blalock  (2003)  used  the  Hvorslev  method,  which 
was  deemed  more  appropriate  (Blalock,  2003;  Entingh,  2002) 

Blalock  (2003)  used  MODPATH  to  determine  a  cumulative  residence  time 
distribution  function.  Fitting  a  polynomial  equation  to  this  distribution,  then  taking  the 
derivative,  gave  the  residence  time  distribution  function  (RTDF).  The  mean  residence 
time,  t,  was  then  calculated  by: 


T  =  /0°°/(0  x  tdt  (8) 

where  f(t)  is  the  RTDF.  The  mean  residence  time,  x,  was  then  compared  to  the 
theoretical  mean  residence  time,  tbar=Vn/Q,  to  determine  cell  volume  utilization  (Blalock, 
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2003).  Blalock  found  that  the  RTDF  was  bimodal,  with  one  peak  at  about  1.25  days  and 


another,  much  smaller  peak  at  around  3.0  days  (Figure  9). 


Figure  9:  RTDF  Calculated  by  Blalock  (2003) 


Blalock  hypothesized  that  this  pattern  could  be  due  to  flow  leakage  or  heterogeneous 
hydraulic  conductivities.  The  mean  hydraulic  residence  time  (x)  of  the  RTDF  was 
calculated  to  be  1.6  days,  which  is  shorter  than  the  calculated  theoretical  mean  residence 
time  (tbar)  of  2.16  days.  Blalock  concluded  that  this  indicated  hydraulic  short-circuiting, 
with  areas  of  stagnant  flow  contributing  to  the  long  tailing  of  the  RTDF. 
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Corbin  (2008)  is  the  most  recent  study  of  the  hydraulics  of  the  uVFTW  at 
WPAFB.  Her  work  was  essentially  a  replication  of  and  follow-up  to  the  work  of  Entingh 
(2002)  and  Blalock  (2003).  Field  sampling  was  conducted  to  determine  values  of 
potentiometric  head  and  hydraulic  conductivity,  and  these  data  were  then  used  in  a 
computer  model  to  provide  insight  into  the  hydraulic  behavior  of  the  uVFTW  and  to 
compare  current  behavior  to  prior  work  (Corbin,  2008). 

Corbin  (2008)  represented  the  wetland  in  GMS  using  a  three-dimensional  grid  of 
fixed  cell  length  and  width,  with  cell  depth  varying  by  layer  (Corbin,  2008).  The  wetland 
dimensions  of  18  m  by  36  m  were  represented  by  a  grid  of  cells  Vi  m  square,  resulting  in 
a  36  x  72  cell  matrix.  The  overall  1.7  m  depth  of  the  wetland  cell  was  divided  into  nine 
layers.  Having  nine  layers  allowed  the  1:1  slope  of  the  walls  of  the  cell  to  be  represented, 
as  deeper  layers  are  smaller  than  shallower  ones.  Table  7  summarizes  the  grid 
construction.  Data  from  Amon  et  al.  (2007)  are  included  for  reference. 
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Table  7:  Grid  Construction,  Corbin  (200S 

0 

Description 

Thickness 

(m) 

Bottom 

Top 

Width  x  Length 
(#  cells  x  #  cells) 

From  Amon  et  al. 
(2007) 

Open  water 
on  surface 

0.10 

1.60 

1.70 

36x72 

Soil/water 

interface 

0.08 

1.52 

1.60 

36x72 

Hydric  soil 
layer  1 

0.19 

1.33 

1.52 

36x72 

0.38  m  hydric  soil 
(layer  A) 

Hydric  soil 
layer  2 

0.19 

1.14 

1.33 

34x70 

Hydric  soil 
layer  3 

0.19 

0.95 

1.14 

34x70 

0.38  m  hydric  soil 
(layer  B) 

Hydric  soil 
layer  4 

0.19 

0.76 

0.95 

34x70 

Hydric  soil 
layer  5 

0.19 

0.57 

0.76 

32x68 

0.38  m  hydric  soil 
(layer  C,  10% 
wood  chips) 

Hydric  soil 
layer  6 

0.19 

0.38 

0.57 

32x68 

Gravel  layer 

0.23 

0.15 

0.38 

30x66 

0.23  m  limestone 
gravel 

(4cm  diameter) 

Datum 

0.00 

0.00 

0.15 

n/a 

0.15  m  sand 

In  the  flow  model,  the  top  layer  was  represented  as  a  constant  head  boundary. 
Corbin  (2008)  had  attempted  to  represent  the  top  layer  as  a  general  head  boundary,  but 
determined  during  model  calibration  that  this  resulted  in  unrealistic  head  contours  a  full 
meter  above  the  wetland  surface. 

The  influent  piping  was  modeled  as  a  series  of  injection  wells  that  each 
introduced  0.192  m3  d1.  The  wells  are  located  in  cells  15-65,  rows  9,  19,  and  29  (Figure 
10). 
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Figure  10:  Influent  Piping  Modeled  as  153  Injection  Wells  in  3  Rows 


Source:  Corbin  (2008) 


3  1 

This  gives  a  total  of  153  injection  wells;  and  a  total  flow  into  the  wetland  of  0.192  nr  d' 
x  153  =  29.3  m  d'  ,  or  5.38  gpm.  This  total  modeled  flow  closely  matches  the  measured 
value  of  5.5  gpm  (Corbin,  2008),  but  is  well  short  of  the  maximum  inflow  rate  of  54.7  m 
d  1  (10.0  gpm)  that  the  pump  is  capable  of  supplying  to  the  wetland  cell  (Amon  et  al., 
2007). 

Corbin  (2008)  used  a  value  of  0.5  for  soil  porosity  for  all  soil  layers,  citing  Amon 
et  al.  (2007).  Amon  et  al.  (2007)  states  that,  based  on  research  conducted  by  others  on 
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the  type  of  soil  used  in  the  treatment  wetland,  the  soils  “might  release  50%  of  their 
volume  as  water.”  Entingh  (2002)  also  estimated  porosity  at  0.53  (Entingh,  2002). 

Hydraulic  conductivity  in  the  horizontal  direction  was  assumed  to  be  isotropic 
(Corbin,  2008).  Corbin  (2008)  assumed  that  since  the  same  type  of  soil  was  used  in  each 
wetland  cell  layer,  that  there  would  be  no  horizontal  anisotropy.  Furthermore,  the 
equations  used  to  calculate  hydraulic  conductivity  (Hvorslev,  1951)  assume  homogeneity 
and  isotropy. 

The  field  measurements  to  determine  hydraulic  conductivity  were  done  in  a 
similar  manner  to  prior  work,  but  with  a  lower  pressure  flush  of  the  well.  This  was  done 
in  an  attempt  to  avoid  the  problems  of  Entingh  (2002)  who  found  that  hydraulic 
conductivity  changed  by  1  or  2  orders  of  magnitude  after  the  wells  were  developed. 
Corbin  (2008)  found  that  with  the  lower  pressure  development,  hydraulic  conductivity  in 
two  wells  increased  by  an  order  of  magnitude,  while  in  another  well  it  decreased  by  an 
order  of  magnitude,  allowing  her  to  conclude  that  developing  the  wells  produced  more 
accurate  values  for  hydraulic  conductivity  (Corbin,  2008). 

Values  for  horizontal  conductivity  were  derived  from  slug  test  results  using  the 
Hvorslev  (1951)  method  and  then  imported  into  the  model  (Corbin,  2008).  Values  for 
vertical  conductivity  (Kv)  were  not  directly  measured,  but  were  assumed  related  to 
horizontal  conductivity  (Kh)  by  the  ratio  Kh  /  Kv  =  1.5  for  all  soil  layers  (1.0  was  used  for 
the  top  open  water  layer).  Corbin  (2008)  found  that  the  model  was  not  sensitive  to  the  Kh 
/  Kv  ratio,  as  the  model  returned  equivalent  outputs  for  ratio  values  from  1.0  to  3.0. 

Model  calibration  was  conducted  by  varying  the  boundary  head  values  while 
holding  other  parameters  constant  until  the  sum  of  the  residuals  between  calculated  and 
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actual  head  values  was  minimized  (Corbin,  2008).  Corbin  (2008)  used  boundary  head 
settings  of  1.7  m  at  the  surface  of  the  water  and  2.9  m  at  the  bottom  of  the  model, 
resulting  in  a  sum  of  squared  weighted  residuals  (SSWR)  of  20.5  m2.  She  also  found  that 
a  SSWR  of  12.0  m  could  be  achieved,  but  the  model  was  then  producing  unrealistic 
velocities  and  water  discharge. 

Corbin  (2008)  found  that  the  uVFTW  was  generally  producing  upward  flow,  with 
piezometers  in  the  bottom  layer  generally  having  the  highest  head  readings  and 
piezometers  in  the  top  layers  having  the  lowest  head  values  (Corbin,  2008).  However, 
the  range  of  head  values  within  each  layer  was  significant,  and  there  was  no  pattern  to  the 
distribution  of  head  values  within  each  layer.  Comparing  these  data  to  the  Entingh 
(2002)  study,  she  concluded  that  the  overall  variations  in  head  measurements  had  not 
changed  significantly  in  the  five  years  between  the  two  studies  (Corbin,  2008). 

Flow  characteristics  were  similar  to  those  found  by  Entingh  (2002).  The  water 
moved  horizontally  in  the  gravel  layer  to  the  region  of  highest  hydraulic  conductivity, 
then  vertically  to  the  surface  (Corbin,  2008).  In  contrast  to  Entingh,  who  found  that  most 
of  the  water  short-circuited  vertically  at  the  north  side  of  the  cell,  Corbin  found  that  the 
water  was  moving  to  the  south  side  of  the  cell,  then  vertically  to  the  surface  (Figure  11). 
She  attributed  this  to  a  change  in  the  actual  behavior  of  the  wetland,  as  opposed  to  an 
error  by  either  herself  or  Entingh  (Corbin,  2008). 
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Figure  11:  MODPATH  Water  Path-Lines  from  Bottom  Layer  (Red 
Dots)  to  Top  Layer  (Black  Dots) 


Source:  Corbin  (2008) 


MODPATH  calculations  were  done  to  determine  RTDFs.  Corbin  found  that  the 
water  spent  an  average  of  0.47  days  travelling  through  the  bottom  layer,  0.48  days  in  the 
middle  layer,  and  0.43  days  in  the  top  layer,  for  a  total  of  1.38  days  (Corbin,  2008).  For 
areas  of  the  wetland  where  flow  was  essentially  stagnant,  particles  of  water  took 
>400,000  days  to  travel  from  the  bottom  to  the  top.  Considering  the  theoretical  residence 
time  (tbar)  of  8.75  days,  water  appeared  to  be  flowing  through  less  than  half  of  the  volume 
of  the  cell  (Corbin,  2008).  Table  8  summarizes  RTDF  findings  of  Entingh  (2002), 
Blalock  (2003),  and  Corbin  (2008). 
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Table  8:  Summary  of  Hydraulic  Residence  Time  Calculations 


Mean  hydraulic 
residence  time 
<j) 

Theoretical  hydraulic 
residence  time 

(tbar) 

Entingh  (2002) 

Three  days 

Not  calculated 

Blalock  (2003) 

1.25  days 

2.16  days 

Corbin  (2008) 

1.38  days 

8.75  days 

To  determine  the  theoretical  hydraulic  residence  time,  tbar,  Blalock  (2003)  and 
Corbin  (2008)  used  the  same  equation: 


Blalock  (2003)  reasoned  that  since  2/3  of  the  inflow  was  assumed  to  be  leaking  out  the 
bottom  of  the  wetland  and  only  1/3  was  actually  transiting  the  wetland,  V  should  also  be 
reduced  by  2/3.  Using  a  porosity  (n)  of  0.27,  volume  (V)  of  651  m3  x  1/3  =  217  m3,  and 
flow  (Q)  of  81.8  m3  d'1  x  1/3  =  27.3  mJ  d'1,  he  arrived  at  tbar  =  2.16  days  (Blalock,  2003). 
Corbin  (2008)  used  porosity  of  0.5,  volume  of  523  m3,  and  flow  of  30.0  m3  d"1  and 
arrived  at  tbar  =  8.75  days.  Clearly,  the  reasoning  of  Blalock  (2003)  is  flawed  and  the 
calculation  of  Corbin  (2008)  is  more  accurate. 

Although  there  is  significant  variation  in  hydraulic  residence  time  calculations 
between  the  three  studies,  all  three  came  to  the  conclusion  that  hydraulic  short-circuiting 
is  a  problem  for  the  uVFTW  at  WPAFB.  In  chapter  3,  a  baseline  model  will  be 
constructed  similar  to  Corbin  (2008).  Various  proposed  engineering  solutions  to  the 
problem  of  hydraulic  short-circuiting  will  be  incorporated  into  the  model  to  determine 
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their  effectiveness  at  improving  performance.  The  ability  to  evaluate  ways  of  improving 
the  performance  of  complex  systems  such  as  uVFTWs  highlights  the  usefulness  of  the 
modeling  approach. 
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III.  Methodology 


3.1  Overview 

The  purpose  of  this  chapter  is  to  explain  the  steps  used  to  evaluate  the  relationship 
between  uVFTW  cell  characteristics  (hydraulics,  layer  depths,  degradation  rates,  etc.)  and 
the  resulting  treatment  efficiency  of  the  cell.  The  first  step  was  to  develop  a  model  that 
simulates  uVFTW  operation.  Model  development  began  by  recreating  the  flow  model  of 
Corbin  (2008)  using  the  MODFLOW  and  MODPATH  modules  of  the  Groundwater 
Modeling  System  (GMS).  The  flow  model  was  then  coupled  with  a  transport  code  that 
simulates  first-order  chlorinated  hydrocarbon  degradation  using  first-order  degradation 
rate  constants  from  the  literature.  After  development  of  the  model,  the  model  was  used  to 
quantify  the  sensitivity  of  simulated  treatment  efficiency  to  varying  design  and 
environmental  parameter  values.  Finally,  model  simulations  were  run  to  determine  the 
effectiveness  of  proposed  engineering  improvements  at  improving  treatment  efficiency. 


3.2  Model  Development 

Groundwater  Modeling  Software  (GMS)  version  7.1.9  was  the  software  package 
used  for  this  research  effort.  It  is  the  current  version  of  the  same  software  used  in  Corbin 
(2008).  GMS  is  a  user-friendly  interface  to  MODFLOW-2005,  a  commonly  used 
program  used  to  model  groundwater  flow  in  three  dimensions.  MODFLOW  is  based 
upon  a  finite-differences  solution  to  the  following  equation  for  groundwater  flow 
(Harbaugh,  2005): 
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(17) 


d 

dx 


dh\  d  /  dh\ 
dy)  +  dz \zz  dz) 


+  W  = 


dh 

dt 


where, 


K  Hydraulic  conductivity  along  the  x,  y,  and  z  axes  [L  T  '] 

h  Potentiometric  head  [L] 

W  Source  (W>0)  or  sink  (W<0),  volume  water  per  unit  volume  per 

time  [T  *] 

Ss  Specific  storage  of  the  media  [L  1  ] 

t  time  [T] 


Since  it  is  generally  not  possible  to  solve  the  equation  for  flow  analytically, 
MODFLOW  uses  a  finite-differences  method  which  is  based  upon  a  finite  set  of  discrete 
points  and  changes  in  head  values  between  these  points  (Harbaugh,  2005).  These  points 
represent  the  center  of  cells,  and  the  model  space  is  thus  represented  as  a  matrix  of  cells. 
The  equation  for  groundwater  flow,  expressed  in  finite  difference  form,  is  the  balance  of 
flows  into  and  out  of  the  cell  set  equal  to  the  change  in  storage  of  the  cell: 


(18) 


where, 


Qi  Flow  rate  into  the  cell  [L3  T  1  ] 
SS  Specific  storage  [L'1] 

AV  Volume  of  the  cell  [L3] 

Ah  Change  in  head  over  time  At 
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MODFLOW  simulates  steady  state  conditions  by  setting  SS=0.  The  equation  for  flow 
then  becomes: 


(19) 


which  simply  states  that  the  sum  of  flows  into  and  out  of  any  cell  must  equal  zero 
(Harbaugh,  2005).  Solving  a  set  of  these  algebraic  equations  for  all  cells  in  the  matrix 
yields  a  numerical  solution  to  the  equation  of  flow.  Using  the  finite  differences  method, 
flow  into  and  out  of  cells  is  calculated  according  to  Darcy’s  law: 


Q  i,j -1/2, k 


KRij-i/z.khCiAv, 


(hi,j-l,k  h'ij.k') 


A  r. 


j~  1/2 


(20) 


where, 


Qi.j— 1/2, k 
KRi,j-i/2,k 

ActAvk 
At j -i/2 


head  at  node  (center  of  cell)  i,j,k  [L] 

flow  rate  through  the  face  between  cells  i,j,k  and  i,j-l,k  [L3  T"1] 
hydraulic  conductivity  along  the  row  between  nodes  i,j,k  and  i,j- 
l,k  [LT1] 

area  of  the  cell  faces  normal  to  the  direction  of  flow  [L'] 
distance  between  nodes  i,j,k  and  i,j-l,k  [L] 


The  formulation  of  the  subscripts  allows  similar  expressions  to  characterize  the  flow 
through  all  six  faces  of  each  cell. 

The  design  of  the  uVFTW  at  WPAFB  lends  itself  to  the  grid  approach  to  model 
design  in  MODFLOW.  The  model  was  designed  as  a  grid  of  36  cells  by  72  cells,  each 
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0.5m  square,  for  a  total  grid  dimension  of  18m  x  36m.  The  1.7m  depth  of  the  wetland 
was  structured  as  10  layers,  with  properties  outlined  in  Table  9. 


Table  9:  Grid  Construction,  Current  Research 


Description 

Thickness 

(m) 

Bottom 

Top 

Width  x  Length 

Kh 

md1 

Kh/ 

Kv 

Porosity 

Layer  1, 

Const  Hd  boundary 

0.10 

(surface  at 
cell  center) 

1.60 

1.80 

cells:  36  x  72 

meters:  18x36 

Interpo 

lated 

1.0 

0.3 

Layer  2, 

Hydric  soil 

0.20 

1.40 

1.60 

cells:  36  x  72 

meters:  18x36 

Interpo 

lated 

1.0 

0.3 

Layer  3, 

Hydric  soil 

0.20 

1.20 

1.40 

cells:  34  x  70 

meters:  17  x  35 

Interpo 

lated 

1.0 

0.3 

Layer  4, 

Hydric  soil 

0.20 

1.00 

1.20 

cells:  34  x  70 

meters:  17  x  35 

Interpo 

lated 

1.0 

0.3 

Layer  5, 

Hydric  soil  /  gravel 

0.12 

0.90 

1.00 

cells:  32  x  68 

meters:  16  x  34 

Interp 

250 

1.0 

0.3 

Layer  6, 

Hydric  soil 

0.20 

.70 

.90 

cells:  32  x  68 

meters:  16  x  34 

Interpo 

lated 

1.0 

0.3 

Layer  7, 

Hydric  soil 

0.20 

.50 

.70 

cells:  30  x  66 

meters:  15  x  33 

Interpo 

lated 

1.0 

0.3 

Layer  8, 

Hydric  soil 

0.18 

.30 

.50 

cells:  30  x  66 

meters:  15  x  33 

Interpo 

lated 

1.0 

0.3 

Layer  9, 

Gravel/soil  mix 

0.10 

.20 

.30 

cells:  28  x  64 

meters:  14  x  32 

25 

1.0 

0.3 

Layer  10, 

Gravel,  influent 

0.20 

0.00 

0.20 

cells:  28  x  64 

meters:  14  x  32 

250 

2.5e08 

Kv= 

250.0 

0.3 

Datum 

0.00 

0.00 

0.00 

n/a 

Layer  Property  Flow  was  used  as  the  flow  package,  while  the  solver  used  was  the 
Preconditioned  Conjugate  Gradient  Method  (PCG2).  PCG2  was  selected  in  favor  of 
Strongly  Implicit  Procedure  (SIP1)  due  to  its  superiority  at  handling  dry  cells  (Harbaugh, 
2005).  For  simplicity,  precipitation  and  evapotranspiration  were  turned  off. 
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The  constant  head  boundary  in  layer  1  was  used  to  simulate  the  effluent  weir 
removing  water  from  the  wetland.  Head  was  set  to  1.75m  in  layer  1,  which  is  0.05m 
above  the  center  of  the  cell  and  represents  that  depth  of  standing  water  on  the  surface. 
Water  entering  layer  1  from  below  would  have  the  effect  of  increasing  head  in  the  layer. 
MODFLOW  simply  removes  this  amount  of  water  from  the  wetland  model  entirely  in 
order  to  maintain  head  at  the  specified  level.  The  drawback  to  this  method  is  there  is  no 
central  point  from  which  to  measure  effluent  contaminant  concentration.  Effluent 
concentration  will  be  estimated  in  this  research  effort  by  averaging  concentration  data 
from  the  layer  1  cells. 

Two  other  methods  of  removing  water  from  the  wetland  were  considered  and 
rejected.  The  first  was  to  model  the  effluent  weir  as  an  extraction  well  in  layer  1.  This 
was  rejected  due  to  the  artificial  effect  an  extraction  well  would  have  on  flow.  The 
second  was  to  model  the  effluent  weir  as  a  drain  cell.  Water  would  flow  out  of  the 
wetland  if  head  was  above  a  specified  level,  but  no  flow  would  occur  if  head  was  below 
the  specified  level.  Although  this  is  conceptually  a  better  way  to  model  the  effluent  weir, 
in  practice  the  drain  cell  proved  incapable  of  removing  a  sufficient  quantity  of  water. 
Since  the  drain  cell  was  located  in  layer  1,  which  was  also  the  constant  head  boundary, 
the  constant  head  boundary  always  removed  more  water  than  the  drain  cells. 

Rather  than  model  the  influent  piping  using  153  injection  wells  arranged  in  three 
rows  as  was  done  in  Corbin  (2008),  the  influent  pipes  were  modeled  as  3  rows  of  highly 
permeable  cells  connected  to  one  cell  designated  as  an  injection  well  (WELl  package). 
The  injection  well  was  set  to  29.3  m  d'  to  reflect  the  total  inflow  (Figure  12). 


52 


Horizontal  hydraulic  conductivity  in  the  influent  rows  was  set  to  2.5e08  m  d"1,  which  is 
six  orders  of  magnitude  higher  than  the  surrounding  gravel  cells  (250  m  d"1). 


Figure  12:  Influent  Piping  Modeled  as  3  rows  of  Highly  Permeable  Cells 
Connected  to  One  Injection  Well 


The  purpose  for  this  change  to  Corbin’s  (2008)  approach  was  to  more  accurately 
model  the  effect  of  heterogenous  permeability  on  what  is,  in  effect,  a  multi-nodal 
injection  well.  If,  for  example,  one  of  the  three  rows  of  influent  piping  was  in  highly 
permeable  soil,  then  more  of  the  total  flow  should  flow  into  that  row,  and  the  other  two 
rows  would  carry  substantially  less  flow.  Modeling  the  influent  as  separate  injection 
wells  does  not  allow  this  phenomenon  to  be  simulated. 

Another  way  to  model  the  influent  piping  would  be  to  use  the  Multi-Nodal  Well 
package  (MNW).  This  package  is  designed  to  simulate  wells  that  are  screened  in  more 


53 


than  one  level.  Unfortunately,  the  current  package  supported  by  GMS  (MNW-l)  does 
not  simulate  horizontal  wells  appropriately,  and  the  new  package  (MNW-2)  is  not  yet 
supported  by  GMS  (Halford  and  Hanson,  2002). 

The  field  measurements  of  head  from  Corbin  (2008)  were  imported  as  an 
Observation  Data  Set.  MODFLOW  reported  the  difference  between  the  calculated  head 
values  and  the  observed  head  values,  within  a  user-specified  standard  deviation  and 
confidence  interval.  This  difference  between  calculated  and  actual  head  values,  the  sum 
of  squared  weighted  residuals  (SSWR),  is  the  objective  function  MODFLOW  uses  to 
measure  how  close  the  calculated  values  are  to  the  measured  values.  The  lower  the 
SSWR,  the  closer  the  calculated  values  are  to  the  measured  values. 

Corbin  (2008)  estimated  values  of  hydraulic  conductivity  from  16  measurements 
that  were  based  on  field  slug  tests  and  then  calculated  using  the  Hvorslev  (1951)  method 
(Table  10). 
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Table  10:  Estimated  Hydraul 

lie  Conductivil 

ty  Values  from  Corbin  (2008 

Label 

X  [m] 

Y  [m] 

Z  [m] 

Kh  [md'1] 

1 

5.2 

2.8 

1.54 

.03915 

2 

4.7 

3.2 

1.16 

.08123 

3 

4.7 

2.8 

0.78 

51.1997 

4 

8.3 

12.6 

1.54 

.0494 

5 

7.9 

13.0 

1.16 

.02025 

6 

7.9 

12.6 

0.78 

3.15456 

7 

17.9 

10.3 

1.53 

.01962 

8 

17.4 

10.7 

1.16 

.43693 

9 

17.4 

10.3 

0.78 

8.19638 

10 

21.1 

7.6 

1.54 

.07378 

11 

20.7 

8.0 

1.16 

.82333 

12 

27.6 

4.8 

1.54 

.46289 

13 

27.0 

5.8 

1.16 

.05619 

14 

27.0 

4.8 

0.78 

51.1997 

15 

30.2 

13.0 

1.16 

.03187 

16 

30.2 

12.6 

0.78 

11.1293 

For  comparison,  Table  1 1  summarizes  typical  ranges  of  hydraulic  conductivity  for 
different  soil  types. 


Table  11:  Typical  Values  of  Hydraulic  Conductivity 


Soil  type 

Kh  [m  d'1],  lowest 

Kh[m  d'1],  highest 

Source 

Silt 

0.0001 

1.7 

(Domenico  and 
Schwartz,  1990) 

Clay 

0.000001 

0.0004 

Sand,  fine 

0.017 

17.3 

Gravel 

25.9 

2592.0 

Soil 

0.004 

0.43 

(Pardue,  2005) 

Considering  that  the  field  measurements  of  Corbin  (2008)  are  all  made  in  the  soil 
layers,  it  can  be  seen  that  the  values  at  the  elevation  z  =  0.78m  are  much  higher  than 
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would  be  expected.  When  these  16  values  are  interpolated  throughout  every  cell  of  the 
3D  grid,  the  effect  of  the  large  values  is  magnified,  depending  on  the  interpolation 
technique  used  (Appendix  B). 

Interpolation  is  a  technique  whereby  calculated  values  are  assigned  to  all  cells 
based  upon  measured  values  of  specific  cells.  Three  of  the  more  commonly  used 
techniques  are  Shepard’s  Method,  Gradient  Plane  Nodal  Functions,  and  Kriging.  All  use 
the  inverse  distance  weighting  method,  which  is  based  on  the  assumption  that  cells 
closest  to  a  measured  value  should  be  influenced  more  by  that  value  than  cells  farther 
away  (Aquaveo,  2010).  They  differ  in  that  Shepard’s  Method  constrains  the  interpolated 
data  set  to  lie  between  the  minimum  and  maximum  of  the  measured  values,  while 
Gradient  Plane  Nodal  technique  produces  values  that  include  minima  and  maxima 
implied  by  the  measured  data.  Gradient  Plane  Nodal  should  be  a  better  technique  than 
Shepard’s  Method,  unless  the  measured  data  include  significant  outliers.  These  outliers 
result  in  minima  and  maxima  that  diverge  even  farther  from  actual  values.  Also, 

Gradient  Plane  Nodal  will  produce  negative  values,  necessitating  truncation  of  the  data  at 
a  minimum  level.  Kriging  differs  from  Shepard’s  Method  and  Gradient  Plane  Nodal  in 
that  the  weights  used  in  inverse  distance  weighting  are  a  function  of  the  measured  values. 
The  variances  of  the  measured  values  are  plotted  against  distance  between  the  points  in 
what  is  known  as  an  experimental  variogram.  The  model  variogram  is  a  best-fit  trend 
line  to  the  experimental  variogram.  The  equation  of  the  model  variogram  is  used  to 
compute  the  interpolation  weights.  Because  the  interpolation  weights  are  mathematically 
determined  by  the  spatial  relationship  of  the  variance  of  the  measured  data,  Kriging  is 
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considered  a  more  accurate  method  of  interpolation  (Aquaveo,  2010).  Kriging  was  the 
only  interpolation  method  used  in  this  research. 

Table  12  illustrates  the  differences  in  interpolated  hydraulic  conductivity  data 
produce  by  the  three  techniques. 


Table  12:  Hydraulic  Conductivity  Estimates  Produced  by  Three  Interpolation 
_ _ _ Techniques _ _ 


[md1] 

Shepard 

Gradient 

Kriging 

Range 

(truncated) 

0.17  to  49.80 

n/a 

-167.5  to  237.9 
0.0001  to  237.9 

-5.66  to  47.0 
0.0001  to  47.0 

Mean 

(truncated) 

7.32 

n/a 

10.03 

22.88 

6.32 

6.35 

Median 

(truncated) 

5.66 

n/a 

1.39 

1.39 

4.83 

4.83 

Std  Dev 

(truncated) 

5.57 

n/a 

50.06 

38.10 

6.37 

6.34 

A  problem  arose  during  model  construction  that  necessitated  significantly 
modifying  the  field  data  used  by  Corbin  (2008).  It  was  found  that  the  observed  head  data 
set  and  the  calculated  hydraulic  conductivity  data  set  could  not  both  be  correct  within  the 
context  of  the  uVFTW  at  WPAFB.  Consider  Darcy’s  Law: 


Q  =  KAVh 


(3) 


where, 


Q  Volumetric  flow  rate  [L3  T'1] 

K  Hydraulic  conductivity  [L  T  ’] 

A  Cross-sectional  area  perpendicular  to  the  direction  of  flow  ( L  ] 

Vh  Head  gradient  [  -  ] 
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For  the  uVFTW  at  WPAFB,  both  Q  and  A  are  constant.  Therefore,  if  hydraulic 
conductivity  is  very  large,  the  head  gradient  has  to  be  very  small.  Using  values  of 
hydraulic  conductivity  interpolated  from  the  measured  values  of  Corbin  (2008),  which 
the  reader  will  recall  were  calculated  based  on  slug  tests  and  the  Hvorslev  (1951)  method, 
and  the  known  values  of  Q  and  A,  it  was  possible  to  estimate  a  head  gradient.  The 
estimated  head  gradient  was  0.006  [m  m’1],  (1.76m  at  the  bottom  and  1.75m  at  the  top  of 
the  wetland  cell)  which  was  much  lower  than  the  measured  gradient.  Since  the  head  and 
hydraulic  conductivity  data  were  incompatible,  and  the  head  data  were  judged  to  be  more 
accurate,  the  observed  head  data  were  retained  and  the  hydraulic  conductivity  data  were 
forced  to  match  by  simply  dividing  the  measured  values  for  conductivity  by  an 
appropriate  constant  before  interpolating.  This  was  done  by  trial  and  error  until  the 
summed  squared  weighted  residuals  (SSWR)  between  computed  heads  and  observed 
heads  was  minimized.  It  was  found  that  when  the  hydraulic  conductivity  values  in  layers 
2  and  4  were  divided  by  15,  and  the  hydraulic  conductivity  values  in  layer  6  were  divided 
by  45,  the  SSWR  was  minimized.  Dividing  the  hydraulic  conductivity  values  by  one 
constant  did  not  produce  a  realistic  head  gradient  and  a  single  constant  chosen  to  match 
the  calculated  and  observed  head  values  in  the  top  two  layers  resulted  in  a  large  error  in 
the  bottom  layer,  and  vice  versa.  This  necessitated  the  use  of  a  second  constant  for  the 
bottom  layer  values  of  hydraulic  conductivity.  Table  13  summarizes  the  profile  of 
interpolated  hydraulic  conductivity  for  the  soil  layers: 
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T able  13:  Interpolated  Hydraulic  Conductivity  Data 


Range 

0.0001  -1.046  md'1 

Mean 

0.147  md  1 

Median 

0.114  m  d  1 

Standard  Deviation 

0.140  md1 

This  had  the  desired  effect  of  producing  a  realistic  head  gradient  between  the  lower  level 
and  the  upper  level  while  also  retaining  the  hydraulic  conductivity  heterogeneity. 

The  interpolated  hydraulic  conductivity  data  set  was  then  populated  using  the 
Layer  Property  Flow  (LPF)  module  within  MODFLOW.  This  had  the  effect  of 
eliminating  the  gravel  layer  (Layer  10),  as  the  original  scatter  data  set  did  not  include 
gravel  layer  measurements.  To  simulate  a  gravel  layer,  a  value  of  250.0  m  d"1  for 
horizontal  hydraulic  conductivity  was  imposed  on  layer  10  after  populating  the  LPF 
module  with  the  interpolated  data  set.  Vertical  conductivity  was  estimated  by  setting  the 
ratio  Kh/Kv  =  1.0  for  all  layers  except  layer  10.  For  layer  10,  setting  Kh/Kv  =1.0  would 
have  resulted  in  overly  high  values  for  Kv  for  the  highly  permeable  cells.  For  the  gravel 
scenario,  Kv  was  set  directly  (not  as  a  ratio)  at  250  m  d'1.  For  scenarios  with  no  gravel 
layer,  Ki/Kv  was  set  to  1.0  for  the  entire  layer  10  except  the  highly  permeable  influent 
cells.  These  were  set  to  Kh/Kv  =  500,000,000  so  that  Kv  =  0.5  m  d"1,  which  produced 
Kh/Kv  ~  1.0. 

Once  the  flow  model  was  developed,  contaminant  transport  was  added  to  account 
for  degradation.  Contaminant  transport  was  incorporated  using  RT3D,  which  is  a 
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submodule  included  within  MT3DMS,  which  is  itself  a  module  of  GMS.  RT3D  is  a 
multi-species  reactive  transport  model  that  is  designed  to  simulate  natural  attenuation  and 
bioremediation  (Clement,  1997).  RT3D  incorporates  seven  different  pre-programmed 
reaction  modules,  to  simulate  the  more  common  contaminants  and  environmental 
conditions.  The  reaction  module  that  will  be  used  in  this  research  effort  is  the 
aerobic/anaerobic  model  for  PCE/TCE  degradation  (Clement,  1997).  This  module 
assumes  first-order  decay  for  all  contaminants,  and  assumes  that  all  decay  happens  in  the 
aqueous  phase  (Clement,  1997).  The  user  must  specify  first-order  anaerobic  degradation 
rate  constants  for  PCE,  TCE,  DCE,  VC,  and  ethene;  as  well  as  first-order  aerobic 
degradation  rate  constants  for  TCE,  DCE,  VC,  and  ethene.  The  reaction  kinetic 
equations  for  aerobic/anaerobic  degradation,  separated  from  the  transport  equations,  are 
summarized  as  follows  (Clement,  1997): 


d[PCE]  _  KP [PCE]  (21) 

dt  RP 

d[TCE]  _  Yt/pKp[PCE]  -  Kt1[TCE ]  -  KT2[TCE]  (22) 

dt  Rt 

d[DCE]  =  Yd/tKt1[TCE ]  -  Kd1[DCE]  -  KD2[PCE]  (23) 

dt  Rd 

d[VC]  =  Yv/dKd1[DCE]-  Kvl[VC\-  KV2[VC ]  (24) 

dt  Ry 

d[ETH]  =  YE/vKvl\VC]  —  Ke1[ETH]~  KE2[ETH ]  (25) 

dt  Re 
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d[Cl\  =  Ylc/PKP[PCE ]  +  Y1c/tKt1[TCE]  +  Y1c/dKd1[DCE]  +  Ylc/vKvl\VC]  (26) 
cLt  Rc 

Y2c/tKT2[TCE ]  +  Y2c/dKD2[DCE]  +  Y2c/VKV2[VC] 

+  Rc 


where, 

[A]  Concentration  of  species  A  [mg  L  '] 

Ki  1st  order  anaerobic  degradation  rate  constant  for  species  A  [T1] 

Ko  1st  order  aerobic  degradation  rate  constant  for  species  A  [  T1 1 

YB/a  Stoichiometric  yield  coefficient  for  the  anaerobic  reductive 

dechlorination  of  A  to  B,  defined  as  the  molecular  weight  of  B 
divided  by  the  molecular  weight  of  A  [-] 

Yl,  Y2  Stoichiometric  yield  coefficient  for  the  anaerobic  (1)  or  aerobic  (2) 

production  of  chloride,  defined  as  the  molecular  weight  of  chloride 
divided  by  the  molecular  weight  of  the  species  which  produces  it  [-] 
Ra  Retardation  coefficient  for  compound  A  [-] 


The  user  can  accept  the  default  Ki  and  K2  values  for  all  layers,  in  which  case  both 
anaerobic  and  aerobic  degradation  would  occur  in  all  layers.  Alternatively,  the  user  can 
specify  different  Ki  and  K2  values  for  each  layer.  For  anaerobic  layers,  K2  would  be  set 
to  zero,  while  for  aerobic  layers,  Ki  would  be  set  to  zero.  It  can  thus  be  seen  from 
equations  (21)  to  (26)  that  anaerobic  decay  is  sequential,  producing  ethane  and  chloride 
as  final  end-products.  Aerobic  decay  produces  chloride  directly  without  intermediate 
products. 

In  Chapter  2  it  was  demonstrated  that  Monod  kinetics  may  be  approximated  as 
first-order  as  long  as  the  concentration  of  the  chlorinated  ethene  substrate  in  the  influent 
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is  at  least  two  orders  of  magnitude  less  than  the  values  for  Ks.  Suarez  and  Rifai  (1999) 

report  the  following  values  for  Ks: 

PCE:  12.00  mg/L 

TCE:  19.00  mg/L 

DCE:  28.00  mg/L 

VC:  23.00  mg/L 

Waldron  (2007)  reported  an  average  influent  concentration  for  PCE  (S  in  equation 

12)  of  46.5  pg/L.  Since  this  is  less  than  the  reported  values  of  Ks  by  approximately  three 
orders  of  magnitude,  1st  order  degradation  is  a  reasonable  approximation.  This 
approximation  also  assumes  that  biomass  density  (B  in  equation  12)  is  a  constant.  In 
other  words,  the  microbial  biomass  population  must  be  at  steady  state.  Given  that  the 
WPALB  VLTW  has  been  in  operation  for  nine  years,  this  seems  reasonable. 

In  order  to  keep  the  model  run  times  within  reason,  sorption  was  assumed  not  to 
occur,  so  R=1.0  in  the  above  equations.  Influent  concentration  was  set  to  1.0  mg  L  1  for 
PCE,  and  0.0  mg  L  1  for  all  other  species.  Model  run  times  were  also  limited  to  30  days  to 
prevent  excessive  run  times.  First  order  anaerobic  degradation  rate  constants  were  0.4  d"1 
for  PCE,  0.014  d'1  for  TCE,  0.004  d"1  for  DCE,  0.004  d"1  for  VC,  and  ,  3.08  d_1  for  ethene. 
First  order  aerobic  degradation  rate  constants  were  0.0  d"1  for  PCE,  0.78  d'1  for  TCE,  0.90 
d’1  for  DCE,  1.85  d'1  for  VC,  and  0.001  d"1  for  ethene.  All  were  taken  from  literature 
reported  values  (Powell,  2010  and  Suarez  and  Rifai,  1999).  It  was  assumed  that  no 
anaerobic  degradation  occurred  in  the  top  two  layers  which  were  assumed  to  be  well 
oxygenated.  It  was  also  assumed  that  aerobic  degradation  occurred  only  in  the  top  two 
layers.  Treatment  efficiency  was  measured  by  percent  reduction  in  total  VOCs  from  the 
influent,  where  concentration  was  comprised  solely  of  0.006  mmol  L  1  of  PCE  (1.0  mg  L-1); 
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to  layer  1,  where  the  average  molar  concentration  across  the  entire  layer  of  PCE,  TCE, 
DCE,  and  VC  constituted  a  proxy  for  effluent  concentration. 

3.3  Sensitivity  Analysis 

Sensitivity  analyses  were  conducted  to  see  how  varying  design  parameters  and 
environmental  parameters  affect  treatment  efficiency.  Specifically,  baseline  values  for 
anisotropy,  layer  thickness,  hydraulic  conductivity,  and  degradation  rate  constants  were 
varied  to  determine  which  parameters  have  the  most  effect  on  overall  treatment 
efficiency.  Because  of  the  significant  model  run  times  for  RT3D,  sensitivity  analyses 
were  conducted  using  MODPATH  residence  times  as  a  proxy  for  treatment  efficiency. 
The  exception  to  this  was  sensitivity  analyses  for  degradation  rates,  which  could  not  be 
run  other  than  with  RT3D. 

•  The  baseline  value  for  vertical  anisotropy  was  Kh/Kv  =  1.0;  and  values  of  1.5,  2.0, 
2.5,  and  3.0  were  run. 

•  Layer  thickness.  The  original  design  called  for  114  cm  of  hydric  soil  based  on  the 
observation  cited  by  Amon  et  al.  (2007)  that  the  root  zone  of  local  fens  extends 
down  to  at  least  100  cm.  It  has  since  been  observed  that  the  root  zone  in  the 
treatment  wetland  extends  into  the  bottom  layer  (personal  communication  with 
Dr.  Shelley).  This  may  be  oxygenating  the  entire  depth  of  the  wetland  cell, 
reducing  the  anaerobic  volume  of  the  bottom  layer.  Increasing  the  depth  of  the 
treatment  wetland  may  increase  the  effectiveness  of  anaerobic  processes. 
Sensitivity  analyses  were  run  with  wetland  depths  of  2.0  m  and  3.0  m  to  test  this. 
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•  Hydraulic  conductivity  was  varied  by  running  every  combination  of  Kh  =  0.2  m  d"1 
in  layers  1-8,  Kh  =  2.0  m  d"1  in  layers  1-8,  Kh  =  25.0  m  d’1  in  layer  10,  Kh  =  250.0 
m  d'1  in  layer  10,  and  Kh  =  1000.0  m  d'1  in  layer  10.  Layer  9  was  held  constant  at 
Kh  =  25.0  m  d'1. 

•  For  degradation  rates,  sensitivity  analyses  were  run  using  mean  values  and  upper- 
range  values.  Minimum  values  were  not  used,  because  these  are  typically  zero 
(Suarez  and  Rifai,  1999).  Table  14  summarizes  values  of  degradation  rate 
constants  that  were  used. 


Table  14:  First  Order  Degradation  Rate  Constants  Used  in  Sensitivity  Analyses 


Parameter 

Variable 

Mean  [d'1] 

Maximum  [d1] 

PCE,  anaerobic 

KP 

0.01 

0.80 

TCE,  anaerobic 

Kti 

0.004 

0.023 

TCE,  aerobic 

Kt2 

0.15 

1.41 

cDCE,  anaerobic 

Kdi 

0.003 

0.006 

cDCE,  aerobic 

Kd2 

0.59 

1.15 

VC,  anaerobic 

Kvi 

0.001 

0.007 

VC,  aerobic 

Kv2 

1.73 

1.96 

3.4  Model  Application 

The  model  was  then  applied  to  estimate  the  effectiveness  of  various  engineering 
solutions  at  improving  overall  treatment  efficiency. 

•  Omission  of  the  gravel  layer.  Corbin  (2008)  showed  most  of  the  water  moving 
laterally  in  the  gravel  layer  to  one  area,  then  up  from  there  (Figure  11).  Omission 
of  the  gravel  layer  allowed  a  test  of  the  hypothesis  that  lateral  hydraulic 
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conductivity  would  be  reduced  and  flow  would  be  biased  upward  from  the  point 
of  injection. 

•  Configuration  of  influent  piping.  With  the  gravel  layer  omitted,  it  may  be 
necessary  to  inject  the  water  via  more  than  3  PVC  pipes  in  order  to  equally 
disperse  the  influent  across  the  length  and  width  of  the  cell  bottom.  Spacing  the 
pipes  lm  apart  would  allow  13  rows  of  influent  piping  to  be  placed  in  the  bottom 
layer.  Model  runs  will  determine  the  impact  of  this  design  change. 

•  Mid-layer  piping.  A  lattice  network  of  3”  perforated  PVC  piping  installed 
midway  between  the  bottom  and  the  top  of  the  cell  may  “capture”  high  hydraulic 
head  in  one  area  and  re-distribute  it  laterally  throughout  the  cell.  The  upward 
flow  may  be  more  evenly  distributed  from  that  point.  Alternatively,  it  may 
exacerbate  hydraulic  short  circuiting  by  providing  a  high  permeability  route  to  the 
area  of  high  hydraulic  conductivity. 

•  Division  of  influent  piping  network  into  zones,  each  one  linked  to  a  surface-level 
shutoff  valve.  Individual  zones  could  be  shut  off  if  it  became  apparent  that 
hydraulic  short-circuiting  was  occurring  in  a  particular  area  of  the  treatment 
wetland. 
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IV.  Results  and  Discussion 


4.1  Overview 

Figure  13  shows  the  hydraulic  conductivity  contours  in  Layer  10  obtained  by 
kriging  the  16  modified  measured  values  of  conductivity.  The  figure  clearly  shows  two 
zones  of  high  conductivity  in  layer  10.  In  cross-sectional  view  A-A,  these  zones  will 
appear  on  the  left  side. 


Figure  13:  Overhead  View  of  Interpolated  Hydraulic  Conductivity  Values, Layer  10 


Using  these  values  for  hydraulic  conductivity,  MODFLOW  calculated  heads  that 
ranged  from  1.75m  at  the  surface  (the  constant  head  boundary)  and  2.40  m  at  depth,  and 
matched  the  observed  head  values  with  a  SSWR  of  221  m“.  An  additional  check  on  the 
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accuracy  of  the  model  was  to  ensure  that  the  flow  budget  was  accurate.  Specifically,  the 


percentage  difference  between  inflows  and  outflows  should  be  less  than  1%.  Also,  total 
inflow  should  equal  both  the  specified  flow  introduced  at  the  injection  well  and  the 
calculated  outflow  at  the  constant  head  boundary  (Table  15). 


Table  15: 

Baseline  Model  Flow  Budget 

In  flows  [mJ  d'1] 

Out  flows  [mJ  d'1] 

Constant  Heads 

0.0 

29.3 

Wells 

29.3 

0.0 

Totals: 

29.3 

29.3 

Percent  Difference: 

o 

o 

31% 

2  2 

The  SSWR  of  221  m  is  significantly  greater  than  the  20.5  nr  achieved  by  Corbin 
(2008),  but  this  is  attributable  to  a  statistical  discrepancy  between  the  two  research 
efforts.  The  data  set  containing  the  observed  head  values  requires  the  user  to  specify 
either  the  standard  deviation  of  the  values  and  the  confidence  level,  or  the  user  must 
directly  specify  the  confidence  interval.  The  computer  then  uses  this  information  to 
decide  how  close  the  calculated  values  are  to  the  observed  values.  It  was  discovered 
during  model  calibration  that  the  standard  deviation  used  in  Corbin  (2008)  was  almost 
twice  the  actual  standard  deviation  as  calculated  directly  from  the  observed  head  values 
(0.51  vs.  0.29).  Furthermore,  closer  inspection  of  the  data  revealed  two  distributions 
(Figure  14).  The  lower  distribution  consisted  of  29  of  the  44  measurements,  had  a  mean 
of  1.82  m,  a  standard  deviation  of  0.09  m,  and  looked  roughly  normal.  The  higher 
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distribution  contained  the  remaining  15  measurements,  had  a  mean  of  2.41  m,  a  standard 


deviation  of  0.05  m,  and  appeared  skewed  left. 


Figure  14:  Histogram  of  44  Observed  Head  Measurements 


C 

<D 
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Source:  Corbin  (2008) 


For  this  research  effort,  the  observed  head  values  for  the  approximately  normally 
distributed  data  set  were  assigned  the  standard  deviation  of  those  data.  GMS  then 
calculated  an  interval  of  2xSD,  producing  a  95%  confidence  interval.  For  the  second  data 
set,  the  broader  criteria  of  Chebyshev’s  theorem  were  applied.  For  non-normally 
distributed  data,  Chebyshev’s  theorem  states  that,  for  k>l,  at  least  (1  -  1/k  )  of  the  values 
will  fall  within  k  standard  deviations  away  from  the  mean  (McClave  et  al.,  2008).  For 
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example,  94%  of  the  values  will  lie  within  4  standard  deviations  from  the  mean  (1  -  1/16 
=  0.94).  For  the  non-normally  distributed  data,  a  standard  deviation  of  0.10  was 
artificially  imposed  upon  the  values  in  order  to  produce  an  interval  of  0.20,  which  is 
equivalent  to  4  standard  deviations. 

GMS  was  able  to  match  the  n=29  set  of  observations  (which  were  higher  in 
elevation)  within  a  95%  confidence  interval.  However,  the  lower  elevation  observations 
could  not  be  matched  unless  the  interval  for  these  observations  was  directly  specified  at 
0.5,  approximately  10  times  the  standard  deviation  of  the  data. 

Use  of  RT3D  to  evaluate  contaminant  degradation  was  severely  limited  by  the 
computational  time  required  for  each  model  run.  Minimum  run  times  for  30  model-days 
were  12-14  hours,  with  some  runs  aborted  when  it  became  clear  that  they  would  take  over 
600  hours.  Despite  these  limitations,  useful  mass  and  concentration  data  for  models  that 
simulated  degradation  were  acquired.  By  inspection  of  Figure  15,  it  is  clear  that  the 
model  is  at  steady  state,  which  allows  a  comparison  of  percent  reduction  in  total  VOCs 
versus  mean  hydraulic  residence  time.  Comparing  percent  reduction  in  total  VOCs  for 
models  of  proposed  design  changes  with  the  baseline  model  allowed  an  evaluation  of  the 
potential  benefit  of  the  proposed  design  change. 
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Figure  15:  PCE  Concentration  Baseline  Model 


PCE  Concentration  by  Layer  [mg  L_1] 


Initial  RT3D  runs  were  conducted  with  mean  values  for  all  degradation  rate 
constants  (Table  14).  It  was  determined,  based  on  high  concentrations  of  total  VOCs  in 
layer  1  and  the  lack  of  correlation  between  percent  removal  and  mean  hydraulic  residence 
time,  that  mean  values  were  too  low  for  the  range  of  mean  hydraulic  residence  times  in 
the  models.  Consequently,  all  RT3D  models  were  re-run  using  degradation  rate  constants 
that  were  an  arithmetic  average  of  the  mean  and  the  maximum  values  (Table  14). 

Sensitivity  analysis  was  conducted  by  varying  parameters  as  outlined  in  chapter  3, 
and  then  comparing  how  either  mean  hydraulic  residence  time  or  contaminant  mass 
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destroyed  was  affected.  Specifically,  the  percent  change  in  the  output  parameter  (the 
hydraulic  residence  time  or  mass  destroyed)  was  divided  by  the  percent  change  in  the 
input  parameter.  The  absolute  value  of  the  resultant  number  indicates  insensitivity  if 
close  to  zero,  and  sensitivity  if  greater  than  1  (Table  16). 


Table  16:  Sensitivity  Analyses  Results 


Input  parameter 

%  change  in  mean  HRT  / 

%  change  in  input  parameter 

Anisotropy 

0.004  to  0.062 

Layer  thickness 

0.978  to  1.094 

Hydraulic  conductivity 

0.029  to  0.033 

PCE  rate  constant 

0.148 

TCE  rate  constant 

0.413 

DCE  rate  constant 

0.006 

VC  rate  constant 

1.379 

The  model  is  insensitive  to  changes  in  anisotropy  and  hydraulic  conductivity,  but 
sensitive  to  changes  in  layer  thickness.  The  model  is  insensitive  to  degradation  rate 
constants,  with  the  exception  of  the  rate  constant  for  VC. 

4.3  Model  Outcomes  of  Proposed  Design  Changes 

The  baseline  model,  while  based  upon  the  uVFTW  at  WPAFB,  is  not  intended  to 
quantitatively  simulate  its  hydraulic  behavior  or  performance  as  a  treatment  reactor.  A 
key  aspect  of  this  study  is  to  demonstrate  how  modeling  may  be  applied  to  investigate  the 
potential  impact  of  design  changes  on  performance. 
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Of  the  seven  different  models  that  were  evaluated,  the  baseline  was  the  worst  in 


terms  of  hydraulic  performance.  At  6.19  days,  the  mean  hydraulic  residence  time  was  the 
lowest,  which  also  produced  the  lowest  cell  volume  utilization  rate  at  65.2%  (Table  23). 
In  a  MODPATH  release  of  approximately  1,000  particles,  88.9%  of  them  had  hydraulic 
residence  times  less  than  the  nominal  hydraulic  residence  time  of  9.50  days  (calculated 
based  on  flow  of  29.3  m  d'  ,  porosity  of  0.30,  and  wetland  volume  of  928  m  ).  The 
residence  time  distribution  function  (Figure  16)  shows  characteristic  short-circuiting 
(RTDFs  for  the  remaining  models  appear  in  Appendix  A).  Degradation  results  for  the 
Baseline  Model  are  summarized  in  Table  17. 


Figure  16:  RTDF  for  Baseline  Model 


Time  [d] 
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Table  17:  Baseline  Model  Degradation  Results 


PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

These  destruction  efficiencies  will  serve  as  the  basis  for  evaluating  the  effectiveness  of 
subsequent  models. 

Omission  of  the  gravel  in  layer  10  (No  Gravel,  3  Pipes  model)  had  a  slight 
improvement  on  hydraulic  performance.  The  mean  hydraulic  residence  time  increased  to 
6.89  days,  while  cell  volume  utilization  increased  to  72.5%.  Removing  the  gravel,  but 
leaving  in  place  the  three  influent  pipes  seems  to  have  had  the  effect  of  “focusing”  the 
head  on  three  narrow  bands  of  soil  centered  directly  above  each  influent  pipe,  decreasing 
mean  hydraulic  residence  time  (Figure  17).  Omitting  the  gravel  in  layer  10  does  seem  to 
have  reduced  the  impact  of  the  zone  of  high  hydraulic  conductivity,  which  is  centered 
above  the  left-hand  influent  pipe  in  Figure  17.  All  three  of  the  spike-shaped  head  profiles 
demonstrate  a  bias  toward  upward  flow.  It  should  be  noted  that  this  model  was  tested 
only  to  evaluate  the  impact  of  changing  this  one  design  parameter,  and  was  never 
recommended  as  a  design  improvement. 
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Figure  17:  Head  Profiles  from  No  Gravel,  3  Pipes  Model  (Section  View  A-A: 

Looking  Down  Longitudinal  Axis) 


Degradation  data  for  the  No  Gravel,  3  Pipes  model  were  improved  over  baseline. 
Layer  1  concentration  was  lower  than  the  baseline  concentrations  for  all  species,  and  the 
reduction  in  total  VOCs  was  improved  from  81.7%  to  89.6%  (Table  18). 
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Table  18: 


Degradation  Data  for  the  No  Gravel,  3  Pipes  Model 


PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L'1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

No  Gravel  3  Pipes 
[mg  L1] 

0.044 

0.044 

0.002 

0.94E-05 

89.6% 

The  next  model,  No  Gravel  12  Pipes,  differs  from  the  prior  model  in  that  12 
influent  pipes  are  used  instead  of  3.  The  improvement  in  the  hydraulic  performance  was 
dramatic,  with  mean  hydraulic  residence  time  increasing  50%  to  9.40  days,  just  slightly 
below  the  nominal  hydraulic  residence  time  of  9.50  days.  Cell  volume  utilization 
increased  to  99.0%.  Instead  of  the  three  spike-shaped  head  profiles  seen  in  the  No  Gravel 
3  Pipes  model,  a  much  more  evenly  distributed  head  profile  is  demonstrated  (Figure  18). 
The  effects  of  the  zone  of  high  hydraulic  conductivity  can  be  seen  slightly  to  the  left  of 
center  in  the  figure. 
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Figure  18:  Head  Profiles  from  No  Gravel,  12  Pipes  Model  (Section  View  A-A) 


Despite  the  improvement  in  hydraulic  residence  time,  the  reduction  in  total  VOCs 
was  not  significantly  better  than  the  No  Gravel  3  Pipes  model  (compare  Table  19). 


Table  19:  De 

gradation  Data  for 

the  No  Gravel,  12  Pipes  Model 

PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L'1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

No  Gravel  12  Pipes 
[mg  L'1] 

0.028 

0.052 

0.003 

1.32E-05 

90.2% 
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The  second  proposed  design  change  was  to  incorporate  a  layer  of  gravel  or 
perforated  PVC  pipes  midway  between  the  top  and  bottom  of  the  wetland.  Such  a  design 
would  create  a  zone  of  very  high  horizontal  hydraulic  conductivity  and  would  possibly 
serve  to  redistribute  differences  in  head  from  areas  of  high  head  to  low  head.  In  the 
Layer  5  Gravel  model,  horizontal  hydraulic  conductivity  in  layer  5  was  set  to 
250,000,000  m  d  ,  a  value  high  enough  to  reflect  a  layer  of  perforated  PVC  pipes  laid 
side-by-side  in  gravel.  Layer  10  was  gravel,  and  three  influent  pipes  were  used. 

The  improvement  in  hydraulic  performance  was  similar  to  the  No  Gravel  12  Pipes 
model.  Mean  hydraulic  residence  time  increased  to  9.44  days  and  cell  volume  utilization 
increased  to  99.4%.  Some  of  the  improvement  in  hydraulic  residence  time  was  probably 
attributable  to  significant  tailing  in  the  RTDF,  which  appears  heavy  enough  to  pull  the 
mean  to  the  right.  This  is  supported  by  the  observation  that  the  percentage  of 
MODPATH  particles  with  hydraulic  residence  time  less  than  the  nominal  hydraulic 
residence  time  was  greater  in  the  Layer  5  Gravel  model  compared  to  the  No  Gravel  12 
Pipes  model,  even  though  mean  hydraulic  residence  time  was  greater  in  the  Layer  5 
Gravel  model. 

Qualitatively,  introducing  gravel  for  Layer  5  did  produce  a  significant  change  in 
the  head  profile  (Figure  19).  It  was  anticipated  that  the  Layer  5  Gravel  could  produce 
this  change  in  head  profile  while  also  decreasing  the  hydraulic  performance  of  the  model, 
which  appears  to  have  happened  for  this  model.  The  inclusion  of  gravel  in  an  actual 
wetland  would  have  required  different  soil  placement  techniques  than  were  used  in  the 
construction  of  the  uVFTW  at  WPAFB.  This  would  probably  have  produced  different 
distributions  of  hydraulic  conductivity  in  the  soil  on  top  of  the  gravel  layer  compared  to 
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the  soil  below  the  gravel  layer,  and  would  have  produced  different  results  than  those 
predicted  by  the  model,  which  assumed  the  hydraulic  conductivity  of  all  layers,  except 
for  Layer  5,  was  unchanged. 


Figure  19:  Head  Profile  for  Layer  5  Gravel  3  Pipes 


Degradation  data  for  the  Layer  5  Gravel  model  are  listed  in  Table  20. 
Degradation  was  better  for  all  species  except  VC,  and  %  reduction  in  total  VOCs  was 
comparable  to  the  No  Gravel  12  Pipes  model. 
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Table  20:  De; 

gradation  Data  for  t 

the  Layer  5  Gravel  M 

lodel 

PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L'1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

Layer  5  Gravel  3  Pipes 
[mg  L'1] 

0.027 

0.051 

0.003 

1.62E-05 

90.2% 

The  next  model,  Layer  5  Gravel  12  Pipes,  combined  what  seemed  to  work  so  well 
in  the  prior  two  models.  Layer  10  gravel  was  removed,  12  influent  pipes  were  used 
instead  of  3,  and  the  Layer  5  high  hydraulic  conductivity  zone  was  incorporated.  While 
the  improvement  in  hydraulic  performance  was  still  significant,  the  model  did  not 
perform  as  well  as  either  No  Gravel  12  Pipes  or  Layer  5  Gravel.  Mean  hydraulic 
residence  time  was  9.16  days,  and  cell  volume  utilization  was  96.4%.  Notably,  however, 
the  percentage  of  MOD  PATH  particles  with  hydraulic  residence  time  less  than  the 
nominal  hydraulic  residence  time  was  the  second  lowest  of  all  the  models,  indicating  that 
the  relatively  low  mean  hydraulic  residence  time  is  probably  due  to  a  few  MODPATH 
particles  with  very  low  hydraulic  residence  times  rather  than  an  overall  leftward  shift  in 
the  RTDF.  This  is  supported  qualitatively  by  the  appearance  of  the  RTDF  (Appendix 
A).  The  effect  of  the  layer  5  gravel  is  still  observed  in  the  head  profile  (Figure  20). 
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Figure  20:  Head  Profile  for  Layer  5  Gravel  12  Pipes 
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Degradation  data  were  better  than  the  Baseline  model  and  similar  to  the  other  two 
models  examined  so  far.  The  concentration  of  VC  went  down  in  this  model,  whereas  it 
went  up  in  the  prior  model  (compare  Table  21). 
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Table  21:  Degradation  Data  for  the  Layer  5  Gravel  12  Pipes  Model 


PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L'1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

Layer  5  Gravel  12  Pipes 
[mg  L'1] 

0.027 

0.053 

0.003 

1.01E-05 

90.0% 

The  final  proposed  design  change  was  to  incorporate  shutoff  valves  within  each 
row.  If  zones  of  hydraulic  short-circuiting  were  identified  during  system  operation,  the 
influent  pipes  under  those  zones  could  be  shut  off,  redirecting  the  flow  to  the  remaining 
open  influent  pipes.  The  models  that  evaluate  this  proposed  design  change,  1  Pipe 
Shutoff  and  4  Pipes  Shutoff,  were  constructed  so  that  an  entire  row  (oriented  along  the 
longitudinal  axis  of  the  wetland),  would  have  to  be  shut  off.  However,  it  is  easy  to 
imagine  any  number  of  configurations,  similar  to  in-ground  sprinkler  system  designs, 
which  would  allow  small  regions  of  the  wetland  to  be  shut  off. 

Shutting  off  entire  rows  was  a  good  test  design  given  the  nature  of  the  hydraulic 
conductivity  heterogeneity  (Figure  13).  The  hydraulic  conductivity  values  showed  two 
zones  of  high  hydraulic  conductivity  -  one  centered  above  row  30,  and  the  other  centered 
above  rows  24-30.  The  1  Pipe  Shutoff  model  shut  off  the  influent  pipe  in  row  30,  while 
the  4  Pipes  Shutoff  model  shut  off  the  influent  pipes  in  rows  24-30. 

1  Pipe  Shutoff  performed  better  than  baseline,  with  a  mean  hydraulic  residence 
time  of  8.62  days  and  percent  cell  volume  utilization  of  90.7%.  4  Pipes  Shutoff  produced 
a  mean  hydraulic  residence  time  of  17.45  days,  and  a  cell  volume  utilization  rate  of 
183.7%.  This  high  value  however,  was  attributable  to  significant  tailing.  A  few 
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MODPATH  particles  had  residence  times  of  130-277  days.  This  indicates  that  the  high 
mean  hydraulic  residence  time  was  due  to  stagnant  or  cycling  flow.  Despite  this,  the  4 
Pipes  Shutoff  model  was  the  best  performer  in  terms  of  the  percentage  of  MODPATH 
particles  with  hydraulic  residence  time  greater  than  the  nominal  hydraulic  residence  time, 
at  36.7%. 

Degradation  performance  of  the  two  models  was  better  than  the  Baseline  model 
for  all  species  and  for  %  reduction  in  total  VOCs.  However,  the  4  Pipes  Shutoff  model 
performed  slightly  worse  than  all  the  other  designs  as  measured  by  %  reduction  in  total 
VOCs  (Table  22). 


Table  22:  Degradation  Data  for  the  Pipes  Shutoff  Models 


PCE 

TCE 

DCE 

VC 

%  Reduction 
Total  VOCs 

Baseline 
[mg  L'1] 

0.067 

0.086 

0.004 

1.42E-05 

81.7% 

1  Pipe  Shutoff 
[mg  L1] 

0.027 

0.053 

0.003 

1.36E-05 

90.1% 

4  Pipes  Shutoff 
[mg  L'1] 

0.033 

0.058 

0.003 

1.02E-05 

88.8% 

A  summary  of  hydraulic  performance  and  degradation  performance  data  for  all 
models  is  shown  below  (Table  23). 
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Table  23:  Summary  of  Hydraulic  and  Degradation  Data  for  All  Models 


Base¬ 

line 

No 

Gravel, 

3  Pipes 

No 

Gravel, 

12  Pipes 

Layer  5 
Gravel, 

3  Pipes 

Layer  5 
Gravel, 

12  Pipes 

1  Pipe 
Shutoff, 
12  Pipes 

4  Pipes 
Shutoff, 
12  Pipes 

nominal  hydraulic  residence 
time  [d] 

9.50 

tau 

[d] 

6.19 

6.89 

9.40 

9.44 

9.16 

8.62 

17.45 

variance 

[d2] 

34.6 

41.6 

56.5 

103.6 

35.2 

30.0 

514.1 

dimensionless 
variance  [  -  ] 

0.90 

0.88 

0.64 

1.16 

0.42 

0.40 

1.69 

%  actual  HRT  <  nominal 

HRT 

88.9% 

83.3% 

72.2% 

77.8% 

68.3% 

75.3% 

63.3% 

%  volume  utilization 

65.2% 

72.5% 

99.0% 

99.4% 

96.4% 

90.7% 

183.7% 

Average  time  layer  1  -4 
[d] 

2.12 

2.72 

2.98 

3.39 

3.95 

3.01 

4.66 

Average  time  layer  5 
[d] 

0.35 

0.46 

0.47 

1.70 

0.75 

0.47 

0.86 

Average  time  layer  6-9 
[d] 

2.25 

3.10 

3.07 

3.64 

3.97 

3.13 

4.85 

Average  time  layer  10 

[d] 

0.57 

0.71 

0.78 

0.91 

0.96 

0.72 

1.38 

Total 

[d] 

5.28 

7.00 

7.30 

9.64 

9.63 

7.33 

11.75 

%  Reduction  Total  VOCs 

81.7% 

89.6% 

90.2% 

90.3% 

90.0% 

88.8% 

95.6% 

Layer  1  Concentration  PCE 
[mg  L1] 

0.067 

0.044 

0.028 

0.027 

0.027 

0.027 

0.033 

Layer  1  Concentration  TCE 
[mg  L1] 

0.086 

0.044 

0.052 

0.051 

0.053 

0.053 

0.058 

Layer  1  Concentration  DCE 
[mg  L1] 

0.004 

0.002 

0.003 

0.003 

0.003 

0.003 

0.003 

Layer  1  Concentration  VC 
[mg  L'1] 

1.42E- 

05 

0.94E- 

05 

1.32E- 

05 

1.62E- 

05 

1.01E- 

05 

1.36E- 

05 

1.02E- 

05 

Note:  “%  actual  HRT  <  nominal  HRT”  is  the  percentage  of  MODPATH  water  particles  in  a  release  of  a 
total  of  1,000  particles  that  had  a  residence  time  less  than  the  nominal  hydraulic  residence  time.  The 
measure  is  useful  in  determining  if  a  RTDF  with  a  high  tau  due  to  significant  tailing  actually  represents  an 
improvement  over  baseline. 


The  MODPATH  average  time  by  layer  data  show  that  the  water  spent  an  average 
of  1.70  days  in  layer  5  in  the  Layer  5  Gravel  model  (Table  23).  This  is  almost  five  times 
the  0.35  days  spent  in  layer  5  for  the  baseline  model,  and  exceeds  every  other  model. 
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Comparing  Figure  19  and  Figure  17,  it  appears  that  the  layer  5  gravel  is  having  the 
intended  effect  of  capturing  and  redistributing  the  three  spike-shaped  head  profiles 
observed  in  the  No  Gravel  3  Pipes  model.  Furthermore,  the  effect  occurs  well  below 
layer  5. 

The  average  time  spent  in  layer  5  in  the  Layer  5  Gravel  12  Pipes  model  was 
anticipated  to  be  close  to  the  1.70  days  observed  in  the  Layer  5  Gravel  3  Pipes  model,  but 
the  value  was  only  0.75  days.  While  this  is  greater  than  the  values  for  the  other  models, 
with  the  exception  of  the  4  Pipes  Shutoff  model,  it  is  less  than  half  the  value  for  the  Layer 
5  Gravel  model.  This  is  due  to  the  much  more  evenly  distributed  head  profile  observed 
in  the  Layer  5  Gravel  model  with  12  influent  pipes  compared  to  the  3  influent  pipes 
model  (the  orange  contour  line  in  Figure  20  barely  extends  into  layer  8,  while  the  same 
contour  line  in  Figure  19  extends  well  into  layer  8).  The  high  value  for  the  Layer  5 
Gravel  3  Pipes  model  was  caused  by  the  time  required  by  the  layer  5  gravel  to  laterally 
dissipate  this  differential  head,  while  the  relatively  low  value  observed  in  the  Layer  5 
Gravel  12  Pipes  model  is  due  to  the  lower  head  differential. 

All  models  achieved  significant  improvements  over  baseline  with  respect  to 
hydraulic  performance.  Mean  hydraulic  residence  time  increased  significantly  (between 
8.62  to  9.44  days,  compared  to  6.19  days  under  the  Baseline  Model),  and  the  percentage 
of  MODPATH  water  particles  with  hydraulic  residence  time  less  than  the  nominal 
hydraulic  residence  time  also  decreased  (between  68.3%  to  77.8%,  compared  to  88.9% 
under  the  Baseline  Model). 

With  the  exception  of  degradation  of  VC  in  the  Layer  5  Gravel  3  Pipes,  all  models 
performed  better  than  baseline  with  respect  to  treatment  efficiency,  as  measured  by  layer 
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1  contaminant  concentration.  Figure  21  is  a  graph  of  percent  reduction  in  total  VOCs  vs. 


mean  residence  time. 


Figure  21:  Percent  Reduction  in  Total  VOCs  vs.  Tau 
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Although  the  expected  jump  in  %  reduction  in  total  VOCs  is  seen,  almost  the 
entire  increase  in  treatment  efficiency  occurs  over  a  very  small  increase  in  mean 
hydraulic  residence  time,  from  6.19  days  to  6.89  days.  Between  6.89  days  and  17.45 
days,  treatment  efficiency  is  essentially  unchanged. 

To  examine  this,  we  make  use  of  Equation  (27),  which  can  be  used  to  predict  the 
course  of  any  first-order  irreversible  reaction,  given  only  the  residence  time  distribution 
function  and  the  reaction  rate  constant  (Clark,  2009). 
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(27) 


C  f 00 

1  — jr^-  =  1—1  exp(—kt)f(t)  dt 
'-in  “'0 


where, 

1  — 2H£:  Reduction  in  contaminant  concentration  C  [-] 

Qn 

k:  Degradation  rate  constant  [d  *] 

f(t):  Residence  Time  Distribution  Function  (T  1  ] 

Equation  (27)  assumes  only  one  contaminant  with  only  one  degradation  rate 
constant,  which  doesn’t  apply  to  the  current  case,  where  sequential  decay  and  multiple 
degradation  rate  constants  are  modeled.  However,  for  purposes  of  illustration,  Equation 
(27)  was  applied  to  the  Residence  Time  Distribution  Functions  for  both  the  Baseline 
model  (tau  =  6.19  days)  and  the  No  Gravel  12  Pipes  model  (tau  =  9.40  days)  with  four 
degradation  rate  constants  (0.001,  0.01,  0.1,  and  0.8  d'1)  which  roughly  correspond  to  the 
range  of  degradation  rate  constants  for  anaerobic  decay  of  PCE  (Suarez  and  Rifai,  1999). 
Figure  22  displays  the  results.  At  very  low  values  for  k,  there  is  no  appreciable 
difference  between  the  two  models  (0.006%  for  the  Baseline  model,  0.009%  for  the  No 
Gravel  12  Pipes  model).  This  is  because  for  the  values  of  tau  used  in  the  two  models,  the 
degradation  time  scale,  as  approximated  by  the  reciprocal  of  the  rate  constant,  is  too  large 
to  show  any  difference  between  the  models.  At  very  high  values  for  k,  there  is  also  no 
appreciable  difference  between  the  two  models  (94.7%  for  the  Baseline  model,  97.0%  for 
the  No  Gravel  12  Pipes  model).  This  is  because  the  degradation  time  scale  is  much 
smaller  than  tau.  It  is  only  when  the  values  of  tau  (6.19  and  9.40  days)  and  the 
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degradation  time  constant  (10  days)  are  similar  that  there  is  a  noticeable  difference  in 
performance  predicted  by  the  two  models  (42.2%  for  the  Baseline  model,  53.4%  for  the 
No  Gravel  12  Pipes  model). 

The  behavior  of  Figure  21  is  now  understood  in  terms  of  the  relationship  between 
tau,  which  ranges  from  6.19  days  to  17.45  days;  and  the  degradation  constant,  which  can 
be  approximated  by  the  value  used  for  anaerobic  decay  of  PCE  in  the  model,  0.40  d'1. 
The  time  constant  for  this  value  is  approximately  2.5  days,  so  an  increase  in  tau  from 
6.19  to  6.89  is  significant,  while  an  increase  from  9.40  to  17.45  is  not. 


Figure  22:  Percent  Concentration  Reduction  vs.  Degradation  Rate  Constant 
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V.  Conclusions  and  Recommendations 


5.1  Summary 

The  goal  of  this  research  effort  was  to  investigate  how  the  hydraulics  of  uVFTWs 
affects  treatment  efficiency.  Specifically,  the  impact  of  hydraulic  short-circuiting  upon 
treatment  efficiency  was  examined.  Several  engineering  design  changes  were  also 
evaluated  with  respect  to  their  impact  on  the  hydraulic  behavior  of  the  uVFTW  and  their 
impact  on  treatment  efficiency.  This  study  confirms  the  existence  of  hydraulic  short- 
circuiting  determined  by  prior  research  by  Entingh  (2002),  Blalock  (2003),  and  Corbin 
(2008).  Furthermore,  this  study  suggests  that  some  simple  design  changes,  such  as 
omitting  the  bottom  layer  gravel  and  introducing  the  influent  via  12  influent  pipes  instead 
of  3,  would  significantly  improve  the  hydraulic  performance  of  the  uVFTW.  The  impact 
on  treatment  efficiency  of  the  improvement  in  hydraulic  performance  is  significant,  but 
occurs  over  a  very  small  increase  in  mean  hydraulic  residence  time,  from  6.19  days  to 
6.89  days.  Further  increases  in  mean  hydraulic  residence  time,  even  approaching  the 
nominal  hydraulic  residence  time,  result  in  minimal  increases  in  treatment  efficiency. 

The  key  insight  of  this  research  is  that,  during  the  design  or  pilot-scale  stage  of 
uVFTW  development,  studies  must  be  conducted  to  quantify  contaminant  degradation 
kinetics  to  determine  an  appropriate  range  of  hydraulic  residence  time  to  ensure  that  the 
uVFTW  operates  at  an  optimum  level.  It  should  be  noted  that  improving  the  mean 
hydraulic  residence  time  into  the  range  where  further  increases  in  treatment  efficiency  do 
not  appear  is  not  necessarily  undesirable.  However,  when  the  design  changes  necessary 
to  bring  about  such  hydraulic  improvements  are  costly,  they  should  be  avoided. 
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5.2  Study  Strengths,  Weaknesses,  and  Limitations 


The  main  strength  of  this  study  was  its  use  of  state-of-the-art  modeling  techniques 
to  simulate  the  potential  impact  of  alternative  uVFTW  designs  on  performance.  The 
wetland  model  enabled  us  to  evaluate  how  design  changes  resulted  in  changes  in 
hydraulic  performance,  and  how  hydraulic  performance  impacted  overall  treatment 
efficiency.  These  results  should  therefore  be  helpful  in  guiding  designs  of  future 
uVFTWs. 

The  main  limitation  of  this  study  was  that,  as  a  modeling  study,  there  is  no 
assurance  that  the  model  reflects  reality.  The  uVFTW  at  WPAFB  and  the  prior  research 
conducted  into  the  degradation  and  hydraulics  inspired  the  current  research,  but  there  was 
no  quantitative  match  to  either  the  hydraulic  performance  or  the  degradation  profile  of 
the  wetland.  Furthermore,  all  models  have  built-in  assumptions  which  may  lead  to 
invalid  results.  Assumptions  such  as  isotropy  may  need  to  be  examined  more  carefully  in 
conjunction  with  a  careful  analysis  of  the  construction  techniques  used  in  building  the 
wetland. 

5.3  Recommendations  for  Future  Study 

The  proposed  engineering  solutions  suggested  in  this  study  required  only 
imagination  to  develop.  Presumably  there  are  many  other  possible  solutions,  and  these 
could  be  tested  in  future  studies  using  the  same  approach  taken  in  this  study.  One  idea  to 
test  is  to  lay  down  sheets  of  impermeable  geomembrane  within  the  wetland,  with  layers 
of  dirt  separating  them.  With  the  proper  design,  water  should  be  forced  to  flow  in  a 
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serpentine  manner  back  and  forth  as  it  moves  upward  through  the  wetland  (think  waiting 
in  line  at  Disneyland).  This  should  increase  hydraulic  residence  time. 

Studies  could  be  conducted  focusing  on  RT3D  to  model  degradation.  RT3D 
allows  incorporation  of  factors  such  as  substrate  limitations,  rhizosome  oxygenation,  and 
oxidation-reduction  conditions.  This  would  require  significantly  more  computational 
power  than  a  typical  desktop  has,  as  well  as  specific  programming  of  RT3D. 

For  any  research  effort  involving  the  use  of  modeling  software,  the  student  should 
seek  out  training  in  the  use  of  that  software  as  early  as  possible.  If  GMS  is  used,  training 
should  be  sought  out  for  MODFLOW  and  GMS,  including  appropriate  sub-modules  of 
both  MODFLOW  and  GMS. 

5.4  Conclusion 

This  study  has  demonstrated  that  some  proposed  engineering  solutions  to  the 
problem  of  hydraulic  short-circuiting  could  be  effective  at  mitigating  the  problem,  while 
also  improving  treatment  efficiency.  The  conclusion  that  can  be  drawn  from  this  is  that 
some  of  the  proposed  engineering  solutions  should  be  incorporated  into  future  uVFTW 
designs.  Omitting  the  bottom  gravel  layer  and  introducing  the  influent  via  many  tightly 
spaced  pipes  rather  than  a  few  pipes  with  more  space  between  them;  and  incorporating 
shutoff  valves  so  that  small  sections  of  the  influent  piping  can  be  shut  down  if  hydraulic 
short-circuiting  is  later  discovered  in  a  particular  area  of  the  wetland;  are  two  specific 
design  changes  that  are  general  enough  to  be  included  in  future  wetland  designs. 


90 


Given  that  the  proposed  engineering  solutions  are  inexpensive  and  easy  to 
incorporate  into  future  designs,  there  is  value  in  conducting  model  studies  of  the 
effectiveness  of  such  proposed  solutions,  regardless  whether  the  outcome  is  positive  or 
negative.  Also,  degradation  kinetics  need  to  be  quantitatively  understood  in  order  to 
determine  an  optimum  range  for  hydraulic  residence  time,  and  to  ensure  that  resources 
are  not  wasted  in  an  attempt  to  improve  hydraulic  performance  where  no  improvement  in 
degradation  performance  will  be  achieved.  Hopefully,  this  study  will  serve  to  guide 
designs  of  future  uVFTWs. 
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Appendix  A:  Residence  Time  Distribution  Functions 


Figure  23:  RTDF  No  Gravel  3  Pipes  Model 


Figure  24:  RTDF  No  Gravel  12  Pipes  Model 
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Figure  25:  RTDF  Layer  5  Gravel  Model 


RTDF 

Layer  5  Gravel 


rH 

i 

0.350 

0.300 

2 

> 

0.250 

u 

c 

OJ 

0.200 

3 

o- 

0) 

0.150 

i- 

LL. 

£ 

0.100 

&_ 

o 

0.050 

z 

0.000 

120.00% 

100.00% 

80.00%  °<u 

> 

60.00%  ra 

3 

40.00%  | 

u 

20.00% 

0.00% 


Time  [d] 


Figure  26:  RTDF  Layer  5  Gravel  12  Pipes  Model 
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Figure  27:  RTDF  Row  30  Shutoff  Model 
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Figure  28 RTDF  Rows  24-30  Shutoff  Model 
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Appendix  B:  Interpolated  Hydraulic  Conductivity  Contour  Plots 


figure  29:  Interpolated  Hydraulic  Conductivity  Contours  (Kriging) 
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Figure  30:  Interpolated  Hydraulic  Conductivity  Contours  (Inverse  Distance  Weighted) 
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Figure  31:  Interpolated  Hydraulic  Conductivity  Contours  (Shepard) 
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